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Summary 

This report deals with the theory, formulation, and solution of compressible three- 
dimensional boundary-layer equations with applications to general swept subsonic or 
supersonic wings in laminar flow. A number of modifications and new features are 
incorporated, based on an earlier general procedure described in NASA CR 4269, Jan. 
1990. A more efficient algorithm has been employed, and overall improvements have 
been made that result in a user-friendly computer code. An interface routine is presented 
that uses the inviscid Euler solutions as input. Code modifications are implemented for 
application in laminar flow control design applications. Output of solution profiles and 
quantities required in boundary-layer stability analysis is included. Conversion routines 
to compare results with Navier-Stokes profiles are also presented. 

This report is a stand-alone document that provides all the necessary details for 
numerical calculation of three-dimensional swept-wing boundary layers. Examples of 
applications and validation with thin-layer Navier-Stokes solutions are presented. A 
user’s manual is included as an appendix. 
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1. INTRODUCTION 


A renewed interest has developed in the design of wings with extensive lengths of 
laminar flow in the subsonic and supersonic regimes. Design for laminar flow by passive 
or active means is a multiparameter optimization problem that involves such variables as 
surface pressure gradients, leading-edge radius, sweep, suction rates, and free-stream 
conditions. To aid in the design process, a reliable computational procedure is needed 
to predict boundary-layer stability. An important part of this computational prediction is 
the accurate generation of smooth mean-flow profiles. This report addresses the issue 
of generating these profiles. 

Two options are available for mean-flow prediction. The first one is the use of an 
accurate thin-layer Navier-Stokes solver in which particular attention is paid to such 
issues as grid resolution and numerical dissipation. However, for repeated preliminary 
design calculations, the Navier-Stokes solution is expensive. The second option is 
the use of an accurate boundary-layer method coupled with an inviscid Euler solution, 
which is particularly attractive for experimentation with different pressure and suction 
distributions. 

This report deals with the theory, formulation, and solution of compressible three- 
dimensional boundary-layer equations, with specific reference to general swept subsonic 
or supersonic wings in laminar flow. A number of modifications and new features are 
incorporated from an earlier general procedure described in NASA CR 4269, Jan. 1990 
(i.e., ref. 1; see also ref. 2). However, the present report is a stand-alone document 
that provides all of the necessary details for numerical calculation of three-dimensional 
swept-wing boundary layers. 

The modifications to the original procedure provide a more efficient algorithm. The 
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solution scheme has been modified to solve the continuity, energy, and momentum equa- 
tions simultaneously, with iterative update of nonlinear terms. Streamwise and spanwise 
differencing schemes have been modified to ensure that the boundary-layer solution 
is consistent with the boundary-layer-edge boundary conditions. Overall improvements 
have been made that result in a more user-friendly computer code. An interface rou- 
tine has been developed to use the Euler solutions as input. Code modifications have 
been implemented for application in laminar flow control design. The modified code also 
provides the output of profiles and quantities that are required in boundary-layer stability 
analysis. Conversion routines that enable one to compare the boundary-layer solution 
profiles with Navier-Stokes solution profiles are also provided. 

Two applications of the code are presented: a subsonic case and a supersonic 
case. For the subsonic case, validations with the thin-layer Navier-Stokes solutions are 
presented. Detailed comparisons are made of the solution profiles and other boundary- 
layer properties. For the supersonic case, comparisons are given of the present code 
with a conical swept wing boundary-layer code developed by Kaups and Cebeci (ref. 3). 

A user’s manual is included as an appendix to this report. The complete program 
package is archived in the NASA Langley computer system mass storage and can be 
made available per individual request. 
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2. FORMULATION 

2.1 Three-Dimensional Boundary-Layer Coordinates 

We start with the definition of the surface-oriented curvilinear nonorthogonal coordi- 
nates used in the three-dimensional boundary-layer equations. 

Let us assume that the surface of interest is defined in terms of Cartesian coordinates 
(x'*, 1 /*, z 1 *) in dimensional units (see Figure 1(a)). The free-stream quantities are 
(A/oo, T£) with a corresponding free-stream velocity U^. 

The body coordinates are normalized with a reference length L the Cartesian 
components of velocity (it'*, v’\ w'*) are normalized by the reference velocity U^. The 
resulting normalized coordinates are ( x\ y\ z'), and the normalized velocity components 
are ( u\ v\ w'). (See Figure 1 (b).) The normalization of other flow quantities is discussed 
in the next section. 

The boundary-layer coordinates are defined in the two surface directions x and y, 
where x is the predominant streamwise direction and y is the general spanwise direction 
(not necessarily along the constant percent chord direction for a wing). The coordinate z 
is mutually perpendicular to both x and y. The coordinates x and y are surface conforming, 
but do not need to be measured as surface arc lengths. For example, for the case 
of the attachment-line flow on a swept wing, the boundary-layer coordinate y may be 
defined along the attachment line on the surface, but can be expressed as distance 
in the spanwise direction perpendicular to the chord. The coordinates x and y can 
also be normalized to the (0, 1) range in the computational domain. The three basic 
metric quantities {hi, h- 2 , 912 ) defined later in this section characterize the stretching and 
shearing of the physical grid into the computational grid. 

The grid line x = 0 coincides with the starting location of the boundary layer as 
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shown in Figure 1 (c) (in this case, the attachment line of a swept wing). The choice of 
x = 0 is important because the system of equations becomes singular at the boundary- 
layer origination point, and a special set of equations must be solved here to initialize 
the solution. 

The transformation of quantities from the (x', y\ z ') to the (x, y, z) system is deter- 
mined by the set of six partial derivatives (x' x , y' x , z' x , x' y , y' y , z' y ). Note that in boundary- 
layer theory, the planes where z = constant are assumed to be parallel to the surface, 
which means that the surface partial derivatives are sufficient to characterize the trans- 
formation from one system to another. These partial derivatives enter into the three- 
dimensional boundary-layer equations via the three metric quantities (hi, h 2 , y 12 ) de- 
fined as 


= \/ K ) 2 + (2/ x ) 2 + (4) 2 (1) 

h 2 = v / ( 4) 2 + fe ) 2 + ( 4) 2 ( 2 ) 

912 = hi h 2 cos /? = x’ x x' y + y' x y' y + z' x z y (3) 

The metric coefficients hi and h 2 represent the stretching in the two surface coordinate 
directions with reference to the physical grid. The nondimensional surface arc lengths 
in the two directions and $ 2 are defined by the two incremental relations 


ds\ = h\ dx 


(4) 


ds 2 = h 2 dy 


(5) 
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Further, the quantity gnlhih 2 is the cosine of the angle j3 between the two coordinates 
x and y. The additional coefficients, based on hi, h 2 and gn, are defined as 


C13 = yjh\h\ - g\ 2 


( 6 ) 


iO 9 12 J 9 12 t . r 

W4 — + "ly 


- ^> 2 *} 


(7) 


C25 = 775- - 2^12/121} ; = ^1*2^ 1 + 7TP7 1 

C 13 l J 


( 8 ) 


g i 2 


C26 = 1 ^ 12 ^ — h 2 h 2 x — ^j^h 2y 


} 


0) 


CW 


r2 l VlZx 
U 13 


<7l2x — h\h\y 


912 , \ 
hi 11 j 


( 10 ) 


C 35 — ,,2 { /( 3' 2x 2^12 ^lj/} 

C 13 


(11) 


C; 


36 


912 

C 2 

l 13 



+ h 2 x 



( 12 ) 


The transformation of velocities from (V, 1 /, «/) to («, t>, u?) is accomplished by 
the inversion of the system given below. Note that it is the symbol used for the 
nondimensional surface-normal velocity; the symbol w is reserved for later use as 
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transformed normal velocity. 


where 


'4 

4 

ill 

V 


U 


~ u f ' 

y'x 

y'x 

£ 2 
0 


V 

— 

V* 

.4 

Zj. 

l — 


_ w _ 


_ w 1 _ 


01 

= y'x 2 

f 

y 

f f 
y 



(13) 


(14) 


02 


(15) 


03 = x' x y' y - x' y y' t 


(16) 


0 



+ 


02 


+ 


,i,2 

¥3 


(17) 


The metric terms defined above allow the flow to be described in terms of surface-oriented 
quantities. As a result, the flow can be treated as equivalent to that on a developed flat 
surface described by two families of nonorthogonal coordinates. However, the metric 
terms do not represent the effects of surface curvature on the physics of the flow. Note 
that the effect of surface shape on the boundary layer is felt only through the inviscid 
pressure distribution. 

2.2 Three-Dimensional Boundary-Layer Equations 

The three-dimensional laminar, compressible boundary-layer equations in surface- 
oriented curvilinear nonorthogonal coordinates are given below (subscripts indicate par- 
tial derivatives). 

The continuity equation is 


6 



+ (ir~ pv) +(^13 pw)x - 0 


The momentum equation in the x direction is 


— U x -f Uy + W Uz + C'24 U 2 + C2*sUV + C*26 V J (M u z) 


p \M Ux 


h\gi -2 


< 2 P , P 

— “x + „ 2 r ! 


C 2 V7M^ 


The momentum equation in the y direction is 


f—v x + JL V + wvi + CW + C'asww + Cse ^ 2 ) - (/^i) 
\/«-l «2 2 


^2912 


P _ P 

l x ^2 r y 


C? 3 * C? 3 V7M^ 


The energy equation is written in terms of the nondimensional total enthalpy H. When 
perfect gas and constant specific heats are assumed, the dimensional total enthalpy H* 


is defined as 


JL^T* + l -q * 2 

7 — I x 


where 


f 2 = u '* 2 + V 1 * 2 + w'* 2 « u * 2 + V * 2 + 2 u* v* cos f3 


The nondimensional total enthalpy H = H* jH^ is based on the reference value 


* R*i 

TT* __ !_ 

00 7 -1 


7 1 * -L — U* z 

1 oo ' 2 00 


which results in the definition 


H = 


T + 0-5 (7 — 1) g 2 
1 + 0.5(7- 
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with T = T*/T^ and 


q 2 = u 2 + v £ + 2 u v cos 


(25) 


The resulting energy equation is 


P (u H x + v H y + w Hg ) = 


±- Hi + iifi-J-W).] 

Pr 2 V Pr J K )z _ 


(26) 


In the equations above, 5 and w are stretched quantities with a tree-stream Reynolds 
number scaling applied as 


= Zy/Rt 


(27) 


w = wy/We oo (28) 

where 

Re - = l&Ll, / i4 (29) 

Other variables are normalized as p* by p T* by T£,; P* by P^; and by p^. The 


equation of state is then written as 


P = P T 


(30) 


The pressure coefficient can be written as 


C v = 


2(^ — 1) 
7 AT 


2 

(X) 


(31) 


The variation of viscosity with temperature is modeled by the Sutherland law as 




p* T * L5 

T* + T* 


T* = 198.6 °R = 110.33 K 


( 32 ) 


= 2.27xl0-»^i = 1.458 x 10 -6 N “ SeC 


f t 2o R l/2 


m 


2 Kl/2 
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In nondimensional form, the Sutherland law is expressed as 


P = 


r La (i + T r ) 

T + T r ' 


Pr . T 

Pr = — , 1 r ~ 

Poo 


1 r 

T* 

± oo 


(33) 


With the assumption that pressure is constant in the normal direction and that the normal 
derivatives tend to zero at the boundary-layer edge (which is denoted by the subscript 
e), the right-hand sides (RHS) of equations (15) and (16) can be replaced by the edge 
quantities to yield 


( U x -f- Uy + WUZ + C24 V . 2 + C25 uv + C26 ) {P U z)i 

\h 1 «2 / 

= p e (^~Ue,x + J[~ Ue ' y ^ 24M e + ^25 U e V e + C26 ^ 


(34) 


( Y~V X + T~Vy + WV' Z + C34 u 2 + C35UV + C^v 2 ) — (pV' z )-. 

\h 1 <12 J 

= p e (-r—V e x + T~ v e,y + Ci 4 ui + C35 U e V e + C36 j 

\ftl «2 / 


(35) 


The equations are hyperbolic in the stream-surface directions and parabolic in the 
surface-normal direction. The boundary conditions (u e , v e , H e ) are specified at the 
boundary-layer edge. The value of w w and one of the values, T w , H w , or // w , are 
specified at the surface. The solution profiles at x = 0 and y = 0 are required to initiate 
the solution marching procedure. 

2.3 Transformation 

A transformation is required for the z and w variables to handle the singularity of the 
equations at x = 0. The transformation also reduces the boundary-layer growth in the 
computational coordinates. Further, a transformation for w is necessary to express the 
resulting equations in a closed form. 
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The transformation is defined as 


C = 


Up 


7 


p dz 


(36) 


t = X\ TJ 


X 

= y] s> = / 


h t 


(37) 


Although £ = x, the partial derivative ^ is actually ^ \ Vi z, and the partial derivative ^ 
is U,c- A similar distinction exists between ^ and The transformations of the 
gradients of an arbitrary function / between the (x, y, z) and (£, y, Q systems are given 
by 

f /d |1 0 **| \f x | 

(38) 


(39) 


' fi 


i — 

uy- 

O 


7/ 

fv 

= 

0 1/,, 


fy 

L/cJ 


o 

o 

N* 

* 


_h. 

7.1 


"1 0 Cd 


r/d 

fy 

= 

o i 


fv 

[fi. 


t 

o 

0 

1 


L/c J 


Because the (3 x 3) matrices are inverses of each other, the following relations are 
valid: 


dC 

dz 



U P 


Pe Pe 5 1 


(40) 


<K 

dx 


HP 


Uf 


Pe Pe 5 1 


(41) 


_ _ ~ / Uc 

dy . * V P \ Pe Pe Sl 


(42) 


We also define the new variables F and G and a transformed normal velocity variable 
w such that 

F = u/u € ; G = v/v r (43) 
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IV 


(44) 


- 51 ^ L 

= U’ — — + 

u e oz 


si_ F #C sj^G^dC 
hi dx ho u e dy 


where v r is an arbitrary reference velocity for v. The transformed velocity w simplifies 
the continuity equation by the explicit removal of the density term. The transformation 
for w can also be rewritten as 


W — 



(45) 


2.4 Transformed Equations 

The application of the transformations for I and w to equations (18), (34), (35), and 
(26) is straightforward and is described in detail in reference 1. (Note that in the present 
report the notations for coefficients have been simplified.) The results are summarized 
here. The continuity equation reduces to 

u>£ = A i + A 2 F + .4 j G r) A\G (46) 


A\ = 


si 

hi 


a Sl r * 

A > = — -t; — < C'13 — 

Cn<j) t hj 


A 3 = - 


S j Vr 
h-2 U e 


A 4 — 


si | C)3 M 


Cn<f> l hi Ue 


(47) 

(48) 

(49) 

(50) 
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<f> = y/pe He Si TF e 


(51) 


Also note that at the boundary-layer edge the transformed normal velocity becomes 

w e£ = A 2 + A 3 Ge,i] + A 4 G e (52) 

At the boundary-layer edge, the partial derivatives of F, G, <j> and the metric coefficients 
in the x or y direction are the same as the corresponding partial derivatives in the £ or 
Tj direction. 

To set up the fourth-order Pade differencing, we introduce three new variables L, M, 
and /, defined as the partial derivatives with respect to £ of the basic variables F, G, and 
H, respectively. The £ momentum equation reduces to 

(IL - wF) c = Bi (F 2 )^ + D>{FG) V + B, F 2 + B 4 F G + B b G 2 + B,9 (53) 

The new variables introduced are defined as 

i = * M .0 ( ; I- H 0 0 - (*); I = (54, 

The coefficients B , are given as 

Bi = -A x (55) 


B ‘2 = 


(56) 


Bs 


U ( e — A 2 + G'24 -Si 

III u e 


(57) 


B 4 = 


■Sl Vr 
hi u% 


U € } 7 ] 


— A.\ + C25 


S\ V r 

U e 


(58) 


12 



Be — C 26 


(59) 


The coefficient Be is evaluated at the boundary-layer edge from equation (53) as 

Be = — + B> G e ,j + B 3 + B 4 G e + B^G^} (60) 

where w e ^ is given by equation (52). 

The 17 momentum reduces to 

(l M —wG)q — Ci(FG)^ -f Co (G 2 )^ + C 3 F“ -f C\ F G + G5G" + C $6 (61) 


The coefficients C, are given as 


Ci = -Ai 


C 2 = -A 3 


G 3 = 


Cm si u t 


C4 = -A-2 + $1 C35 + V r £ 

hi v T 


S] Vt 51 

C5 — — A4 + C36 + 7 

U e h‘2 Ue 


The coefficient Ce is evaluated at the boundary-layer edge from equation (61) as 

Ce = - { w e,c Ge + Ci G e , ( + C 2 {G;) v + C 3 + G 4 G e + C 5 G 2 } (6 
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The energy equation in terms of the transformed variables becomes 


( 


l 

~Pr 


I 



D\ {F H)^ + D 2 {GH) v + D 3 FH + D 4 GH + D 5 



( 68 ) 


The coefficients D, are given as 

D\ = — A\\ D 2 = —.43; D3 = —A 2 ; D.\ — —A\ (69) 


The term D 5 can be rewritten with equation (25) as 


D 5 


u£ ( l ~ Pr \ 

^ V Pr ) 


i\FL + v 2 r GM + u e v r (FM + GL) | (70) 


2.5 Quasi-Two-Dimensional Equations for Initial Conditions 

The solution of the transformed three-dimensional boundary-layer equations requires 
the specification of profiles of F, G, and H at the £ = 0 and rj = 0 planes as initial 
conditions. In the ideal case, a plane of symmetry flow or a conical flow will exist at 
these planes. As a result, the equations simplify to quasi-two-dimensional form (i.e., 
either £ or 77 derivatives reduce to zero). This simplification permits the solution and 
generation of initial solution planes independent of the full three-dimensional flow region. 

In the case of flow past a wing attached to a fuselage, certain assumptions are 
necessary for simplification of the initial profiles. For a swept wing with a nonsymmetrical 
chordwise wing section at an angle of attack, the attachment-line boundary layer (at £ 
= 0 ) is a true symmetry plane only if we assume that the coefficient C 26 = 0. This 
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assumption is reasonable if the curvature of the attachment line is not significant. A 
more drastic assumption is necessary for the rj = 0 plane. We restrict the calculation 
to a certain region on the wing, where the boundary-layer assumptions are valid. At 
the 77 = 0 boundary, we assume that the rj gradients of the flow variables are equal to 
zero. However, we permit variation of the metric terms as well as the boundary-layer 
edge quantities in both the £ and rj directions. With this assumption, the equations 
again reduce to quasi-two-dimensional form. This assumption is a more general version 
of the conical flow assumption and is termed as the locally infinite, swept-wing (LISW) 
assumption. The profiles at the £ = 0, 77 = 0 location are generated by assuming that 
the attachment line is infinite swept locally as well. 

The equations that are valid for an LISW flow are obtained by setting the coefficients 
A 3 (hence, B 2 , C 2 , and D 2 ) to zero. In the finite-differencing scheme, the ^-direction dif- 
ferencing coefficients are also set to zero. To obtain smooth solutions near this boundary, 
the full three-dimensional solution is relaxed to the LISW solution by incorporating a fac- 
tor u =u(n) to be applied to these coefficients for a few planes adjacent to the rj = 0 
plane. A value of w = 0 gives the LISW solution, whereas a value of u> = l gives the 
full three-dimensional solution. A similar assumption of LISW flow is used at the r? = 
77 max boundary as well. 

The quasi-two-dimensional equations that are valid for an attachment line are more 
involved because the £ momentum equation becomes singular at £ = 0 by substituting 
u = 0. The gradient du/d£ is finite, however. The effect of the strong inviscid flow 
acceleration f e = du e jd( on the attachment-line boundary layer is characterized by taking 
the i derivative of the £ momentum equation. The new definitions for the transformed 
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normal coordinate and normal velocity are 


C = 


fe 

Pe Pe hi 


z 

/ 


p dz 


(71) 


. h\ ( ~ v , 

w = p ' l JFFTeV~h z - 


(72) 


With these special transformations, the transformed equations are in the same form as 
before (equations (46), (53), (61), and (68)) except that the following coefficients are 
changed: 

= o (73) 




(74) 


M — 


h\ v 


1 V r 


h2 fe 


(75) 


Aa = 


C 13 


{«'&) 


(76) 


$ — \/ Pe Pe h] f e 


(77) 


Bi = 2 


(78) 


B, = ^ /«,, - A t + C25 h ‘ Vr 


hlf! 


fe 


(79) 


Ft - r kl V r 
Je 


( 80 ) 
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C 3 = 0 


(81) 


C 4 = -A 2 (82) 

C S = - A, + C« ^ (83) 

Je ^2 Je 

Ci = |( 1 rr)^ ( ' GM, ( (84) 


In the finite-differencing scheme for the attachment-line solution, the ^-direction 
differencing coefficients are also set to zero. 

As before for the case of the LISW equations, the equations that characterize the 
flow at £ = 0 , 7 / = 0 are obtained from the attachment-line equations with the added 
condition A 3 = 0. 

2.6 Equations in Vector Form 

To apply the fourth-order Pade differencing formula to the set of equations, the 
equation set is expressed in vectorial form for convenience. The vector Q is defined 
as 

Q = fw , IL — u;F, IM — wG , l p I — wH , F, G , (85) 

where 

<p - j; ( 86 ) 

The normal derivative Q' = dQ/d ( (note that the prime denotes the partial derivative with 

respect to () can be written from the boundary- layer equations. Before the expression 

for Q' is presented, the following equations are given. 
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The density ratio 9 can be written from equation (24) in terms of the elements of Q as 


6 = E 4 F 2 + E 2 FG + E 3 G 2 + E 4 H 


77 X2 2 TP n \2 a 

F\ — U £ ] t 2 — —l U e v r cos p 

X3 X3 


77 X2 2 r? X 1 

£3 = — — u r ; £4 = — 


X3 


X3 


7 - 1 


xi = 1 + X2 = xi - i; X 3 = ftxi 


X2 <?e 


(87) 


The viscosity ratio / can be written from equation (33) as 

, _ ./p{l+Tr/T e \ 

’ - ve \rrrn) 

The normal derivatives /' and O' can be derived in the form 


( 88 ) 


O' = 2ExFL +E 2 (FM + GL) +2E 3 GM + E 4 I 


(89) 


t = * - 


VO 1 2 V 0 (1 +T, 


r r /T e ) } 


(90) 


The Q\ element is obtained from equation (46). The Q 2 element is obtained by 
substituting for 0 from equation (87) into equation (53). 

Q ' 2 = Bx (F 2 ) f + B 2 (FG) v + B 3 F 2 + B 4 FG + B S G 2 + B 6 H 

(91) 

B 3 = B 3 + Ex B 6 ; B 4 = B,x + E 2 B 6 - B s = B 5 + E 3 B 6 ; B 6 = B 6 E 4 


The element Q 3 is similarly obtained as 


Q' 3 = Cx (FG) n + C 2 (G 2 )^ + C 3 F 2 + C 4 FG + C,G 2 + C & H 


Cz — C 3 + ExCs; C 4 — C 4 + E 2 Cq ; C 5 = C 5 + £ 3 Ce ; (?6 = C^E 4 


(92) 
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The second derivatives V = F", M' = G", and /' = H" required subsequently in 
the vector representation are obtained by differentiation. For example, differentiation of 
Q 2 ' yields 

Q ' 2 = IL ' + LI' — wL — Fw — RHS of equation (91) (93) 

If vf is substituted from equation (46) and rearranged, then 


U = 


1 ■ 
7 . 


B x (F 2 )^ + B 2 (FG) v + B 3 F 2 + Bi FG + B 5 G 2 + B 6 H 


—LI 1 + wL + A\FF je + A^FG^ ; 


(94) 


S3 = S3 + A 2 'i B 4 — S4 + A4 


The expressions for A/ and f are similarly obtained as 


M - 7 


C,(FG), + C 2 (G 2 )„ + C 3 F 2 + C,FG + C 5 G 2 + C,H 

-Ml’ + t vM + A,GF ( + .4 3 GG„]; (95) 


C 4 = C 4 + Ao', C b = C b + A\ 


V = i [Di {FH)^ + D 2 {GH) v + D 3 FH + D 4 GIL + D b 


l y. 


-Il'p + wl + A\HF Jc + AiHGr, ; 


( 96 ) 


D 2 = D b + A-i\ D 4 = D 4 + A 4 


The expression for S 5 can be obtained from equation (70) as 
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D 5 = l [Aj (FV + L 2 ) + A 2 (GM 1 + M 2 ) + A 3 (Fi\F + ML + GL' + LM) 


+ /' 


A] FL + A 2 GM + A 3 (FM + GL) 


Ao = 


U£ 1 - Pr 
HSo Pr 


Ai — Aor^; A 2 = Ao v 2 ; A 3 = Aou e t7 r cos 8 


( 97 ) 


Note that A 3 = A 3 = 0 for the attachment-line equations. The final vector repre- 
sentation is given as 


Q = 


/ w \ 

IL - wF 
IM — wG 
l P I - iv H 
F 

G 

H J 


(98) 


Q' 


A\F{ + A2F + AjGj, + A 4 G \ 

B 1 (F 2 ) f + B 2 (FG\ + B3F 2 + B 4 FG + B 5 G 2 + BqH 
Ci{FG\ + C 2 (G 2 ) v + C 3 F 2 + C 4 FG + C 5 G 2 + C 6 H 
Dx{FH\ + D 2 (GH) v + D 3 FII+D 4 GII + D 5 

L 

M 

I J 


(99) 


Q” = 


( A\L^ + A 2 L + A 3 M V t A 4 M \ 

2 B x {FL\ + B 2 (FM + GL\ + 2 B 3 FL + B 4 (FM + GL) + 2 B S GM + B 6 T 
Ci ( PM + GL\ + 2 C 2 {GM) n + 2 C 3 FL + C 4 (FM + GL) + 2 C,GM + C c / 
D,(FI + HL\ + D 2 (GI + HM) n + D 3 (FI + HL) + D 4 (GI + II M) + D' 5 

L' (equation (94)) 

M' (equation (95)) 

V V (equation (96)) / 

( 100 ) 
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3. DISCRETIZATION 


3.1 Differencing Formulas 

We apply the fourth-order Pade differencing formula in the normal direction. This two- 
point compact scheme is defined in terms of the variable and its two higher derivatives. 
In the present case, if we assume that the indices in the two surface directions (i,j) 
remain constant and that k is the normal direction index, then the discretization at the 
midpoint of k and (k - 1 ) is written as 

Qk — Qk- 1 — —j~ ( Q'k + Q'k- 1) + ~j7T (Qk ~ Qk-l) + ^(^C 5 ) = 0 (101) 


A£ = (k — Cfc— i 


(102) 


The differencing in the surface directions £ and rj are second order (or first order in 
some regions). For example, in the £ direction, if we assume that the indices (j, k) are 


fixed, then 


Q^ — ni Qi T 


( 103 ) 


{aQ} — CI2Q1 - 1 + (13Q1-2 

Similarly, in the rj direction, if we assume that the indices k are fixed, then 

Q n = biQij + 


( 104 ) 

{ bQ } = 62 Qi,j - 1 + b^Qij —2 + h Qt-i,j + hQi-i,j+i 

The value of the coefficients a, and b, are dependent on the location of the point (i,j) in 

the streamwise and the crossflow direction. At present, we restrict the discussion to the 
discretization in the £ direction and will use the short notation given above for the £- and 
//-direction differencing. (See “ Discretization in the (i,j) Directions” for more details.) 


We define a solution vector {S'} = {«?, F, < 7 , H, L, M, I} T . Because the equations 
are nonlinear in {S}, Newton linearization is used to convert the system to a linear matrix 
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inversion problem. If superscript n denotes the current iteration stage, let us define {<55} 
as 

{SS} = S n -S n ~ 1 (105) 

A linear system is now set up to solve for {<55} in terms of the solution at iteration level 
n - l. For example, a term that involves (F n ) 2 is written as 

(F n ) 2 = (F B_1 + SF) 2 « (F" _1 ) 2 + 2F" _1 <5F (106) 

In what follows, the superscript n- 1 is dropped and is taken to imply the known values 
of {5} at iteration n - 1. A few examples of the linearized formulas with this notation 
are given below: 

(F n ) 2 = F 2 + 2F <5F 
F£ = Ft + arfF 

F” = F v + b\8F (107) 

(FGOJ = a\G5F 
F n F£ = FF^ + 8F(a\F + F^) 

3.2 Linearized system at (i,j) 

The system is explicit in £ and rj because of our choice of the finite-differencing 
scheme and is implicit in the surface-normal direction. The linearized system is repre- 
sented at location ( i,j ), which corresponds to the solution at iteration level n as 

“L.'F-] {«*} = {--f} (io8) 

where a\ m and b k lm are elements of the (7 x 7) blocks in the diagonal and superdiagonal 
locations of the linearized block bidiagonal system. The superscript k denotes that the 
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discretization corresponds to the midpoint of k and (k — 1) points. The index l varies 
from 1 to 7, depending on which element of equation (101) is being discretized; m varies 
from 1 to 7, depending on which element of {<55} it multiplies. The (7 x 7) blocks a k m j 
and j6f m are the only nonzero blocks in the system above because we have a two- 
point compact scheme in the k direction. The (7x1) vector {r[‘} corresponds to the 
residual of equation (101), based on the solution {5} at iteration (n — 1). 

For example, examine the discretization of the second element of the system repre- 
sented by equation (101), which can be written as 

(IL - wF) k - (. IL - wF) fe _! 

_ + b 2 (fg) v + b 2 f 2 + b a fg + b 5 g 2 + b 6 h~\ k 

+ [Bi F 2 + B 2 {FG\ + B 2 F 2 + B 4 FG + B b G 2 + B 6 h] 

+ B X {FL\ + B 2 {FM + GL\ + 2 B 3 FL + B^FM + GL ) + 2 B 5 GM + B e l] 

- [25i {FL)z + B 2 (FM + GL) V + 2 B 3 FL + B A (FM + GL) + 2 B 5 GM + i? 6 /] ^ j = 0 

(109) 

where 

Aa = a - C*-i ( 11Q ) 

We can now construct the elements of the blocks af m , b k lm , and {rf } for the case 
l = 2 from above by using the linearization procedure explained earlier. Some examples 
are given below. The coefficient a k j is the coefficient of Swk-i from the second element of 
equation (101), which is discretized at (A- - \), and b k x is the corresponding coefficient 
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«21 - Fk - 1 


hi = ~Fk 

a\ 2 = - =^(2fliaiF + + 2B 3 F + B 4 g) 

^22 = -«»fc - ^(‘2i?iaiF + B,b x G + 25 3 F + B 4 g) 

AC 2 / 

«25 = — /*-i - -pr (25 1G1 F + B 2 6!G + 2£ 3 F + 5 4 G) 

12 V /*-l 

r 2 = negative of LHS of equation ( 109 ) 

Similar expressions can be derived from the 7 elements of the system that are 
represented by equations (98)-(100). These expressions for the 49 elements of a\ m , 
the 49 elements of b k l m , and the 7 elements of { rf } are given in appendix A. 

3.3 Boundary Conditions 

The boundary conditions at the boundary-layer edge (k = he, ( = ( e ) are given as 

F k=ke = 1; Gk-kt = G e = — ; Hk = ke = He (1 12) 

V r 

The boundary conditions at the wall (k = 1, ( = o) are given as 

Fk = i = Gk = i = 0 

H k =\ = Hi or I k= i = I\ ( 1 13 ) 

w k=\ = Ww 

The transformed normal velocity at the wall {io k =\ or w\ or w w ) is specified from equation 
(45) as 

/ * * * \ 

U l = ( ) V’w \/~Rc oo (114) 

\ r oo c oo / 
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where is defined for general three-dimensional flow or attachment-line flow as 

V>„=,p^ or ^ = (115) 

V PeVeUe \ PePeJe 

Suction rate is usually specified in terms of the suction-rate momentum value scaled by 
the free-stream momentum as q s = p\w\! p^U^. 

For an adiabatic wall, h is specified as equal to zero. For a nonadiabatic case, if 
the wall heat flux is known in dimensional units as q* w , then the wall heat flux equation 
can be written as 

r . die) 

q “ ( 7-1 )Pr ,w \dz') w 

If we apply the transformation formulas and note that at the wall 0' w = E±I\ (from equation 
(89)), then we obtain 

7 R* (T^n\ 


• * 
r lw 


(7 - 1 )-p* 


~R&oo ifrqTepe I\ 


(117) 


(118) 


where 0 , is defined for general three-dimensional flow or attachment-line flow as 

, / «e , I Te 

*Pq=\ Or V>q= \ T 

V PePe s 1 V P t ^ ekl 

Application of the heat flux boundary condition involves the specification of h with 
equation (117). A negative value of q* w corresponds to a cold-wall situation (i.e., heat 
flow in the negative z direction). If the wall temperature T w is the known quantity, then 

Hi is specified (from equation (24)) as 

T w 


Hi = H 


W 


(l + Ijl M*,) 


(119) 


The incorporation of the boundary conditions into the linear discretized system results 
in a shift of the rows downward by four. This produces a system that is block tridiagonal. 
The block tridiagonal system is defined by 


;... a k 0 k 7 * 


\8S k 


(. k = 1 , 2 , Ice) 


( 120 ) 
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where 



(121) 


Note that the index k on a k , j3 k , and 7* corresponds to the location of a particular block 
in the ( ke x ke) system of (7x7) blocks. The index k of each of the elements in a 


block refers to the discretization location. The subdiagonal, diagonal, and superdiagonal 


blocks at row location k are obtained as 


k 

nk 

k 

r,k 

fc 

k 

k 

3 11 

a \2 

a 13 

a \4 

a 15 

a 16 

a 17 

n k 

k 

„k 

„k 

k 

k 

k 

2 21 

a 22 

a 23 

a 24 

a 25 

a 26 

a 27 

k 

n k 

k 

k 

k 

k 

k 

2 31 

a 32 

a 33 

a 34 

a 35 

a 36 

«37 

k 

k 


„k 

h 

k 

it 

l 4\ 

^42 

°43 

a 44 

a 45 

a 46 

a 47 

0 

(f 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


b k 

°n 

b k 

°12 

b k 

°13 

a* 

"14 

A* 

a 15 

^16 

f>j 7 

b k 

°21 

h k 

°22 

a 23^ 

b k 

"24 

b k 

°2b 

^26 

b 27 

b k 

°31 

b k 

°32 

^33 

A* 

"34 

hk 

°35 

a* 

"36 

b 37 

b k 

°4\ 

b k 

y 42 

^43 

A* 

"44 

b k 

°45 

A* 

°46 

b 47 

a 51 , 

a k ^ 

a 52 

«53 +1 

"54 

«55 +1 

«56 +1 

a k J l 

a fil 

a Ui 
a 62 , 

«63 +1 

G *+ 1 

"64 

4 +i 

«66 

a*" 1-1 

a 67 

„K + 1 

,T + 1 

„Jfc+l 

J+i 

„fc+l 

E+l 

T+l 

L a 71 

a 72 

a 73 

a 74 

a 75 

a 76 

a 77 

" 0 

0 

0 

0 

0 

0 

0 ' 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

b li l 

Ht 1 

6 53 +1 

b k f 

4+> 

b U l 

6* +1 

°61 

5^+1 

w 62 

6 +1 
°63 

6 

y 64 

b i + 1 

a 65 

6 +1 
a 66 

jl+l 

w 67 

L°71 

J+l 

°72 

4*+i 

°73 

A*+l 

O 74 

J+l 

D 75 

hk+ 1 
°76 

A*+l 
"77 J 


( 122 ) 


(123) 


(124) 
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The boundary conditions at the boundary-layer edge result in a diagonal and RHS 


block as 


P ke \ = 


b k 

°11 

^21 

& 

b k 

°41 

0 

0 

0 


b k 

“12 

b k 

°22 

h k 

“32 

b k 

°42 

1 

0 

0 


^13 

h k 

“l4 

b k 

“15 

b l6 

^17 1 

a 22 b 

b k 

“ 24 

b k 

“25 

& 

*2 7 

& 

b k 

“ 34 

b k 

“35 

^36 

& 

b k 

“43 

b k 

“44 

b k 

“45 


b k 

“47 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 . 


1 

6 ke 


>n 

r k 

r 3 

r 4 


L 


0 



0 



i 

0 

1 


; k = ke 


(125) 


The diagonal block fa and RHS block Si that result from a specified wall heat flux 
condition (adiabatic or otherwise) are as given below (which corresponds to 8I\ = 0): 


l* 1 ] = 


1 

0 

0 

0 

fc+1 


a 


t+1 


«71 


+ 1 


0 

1 

0 

0 

G fc+1 

a 52 

J+l 

fl 62 

J+l 

a 72 


0 

0 

1 

0 

k + 1 
x 53 

J+l 

x 63 

J+l 

x 73 


0 

0 

0 

0 


a 


s+l 


a, 


+1 


l 74 


0 

0 

0 

0 

7^' + 1 

f+l 

J-+i 

*75 


0 

0 

0 

0 

7 fc+ 1 
z 56 
t+1 


qx+1 

“76 


0 

0 

0 

1 

a fc+1 
% 


$ 


G 67 
a *+! 
u 77 



1 — 
o 


0 


0 

1 

"57 

II 

0 

r fc+l 

Ui 

r Ui 

L r 7 -1 


fc=l (126) 


For a specified wall temperature condition, the system becomes 


ip') = 


l 

0 

0 

0 


I 1 


Oc 


+i 


0 

1 

0 

0 

7 fc + 1 

x 52 

fc+1 


0 

0 

1 

0 

,*+1 

3 3 +l 


0 

0 

0 

1 

a k+1 

“54 

alt' 


a 7 


iv. > ;n. 

“72 


71 


*73 


a 


74 


0 

0 

0 

0 

Q fc+i 

°55 

J+l 

“65 

X+l 

a 75 


0 

0 

0 

0 

,*+1 


? 6 h 

f 


: + l 


+ 1 


x 76 


0 ' 


0 ' 

0 


0 

0 


0 

0 

; M = 

0 

fc-j-1 

r fc+l 

a 57 

J+l 


3+1 

a 67 

J+ 1 

a 77 - 


Ui 

L r 7 J 


;k = l (127) 


Appendix A summarizes the information required to construct the block tridiagonal 
system. Implementation of the boundary conditions also requires an update that is based 
on the current solution, which is done with the update of the nonlinear terms. 


Because the Pade formula is a compact scheme based on the solution variables and 
their derivatives at two points that span the local cell center, a stretched grid can be 
employed without degradation of the fourth-order accuracy of the method. A stretching 
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constant k s is defined to exponentially stretch the grid in the £ direction as 


Ck — Cke 



( 128 ) 


3.4 Discretization in the ( /,/ ) Directions 

In the fully three-dimensional region (away from the attachment line and side bound- 
aries, j = 1 or j = nt/lim), the differencing in the surface directions £ and ?? is done to 
second-order accuracy. In accordance with the parabolic nature of the equations, the £ 
derivative is obtained by the three-point upwind-differenced formula 

(f()i ,j ~ Q 1 fi,j + a 2fi~l,j + 

ai = (A£“ - A£?_,)/A; ci 2 = -A£“/A; a 3 = A£?_j/A (129) 

^ = + A£,_i); A£, = — £,- 1 ^; A£ t _i = £,_i j — £,_ 2 j 

When i = 2, the first-order formula with just two points is used, which results in the 
coefficients 

a 1 = !/A6 = 1/(6 - 6); «2 = -ay, rt 3 = 0 (130) 

A function w,- is used to blend the first-order and second-order formulas in a small number 
of marching steps. The short notation {aF} used in equation (103) can thus be expanded 
in terms of the coefficients given above. At £ = 0, the attachment-line equations do not 
contain any djd( terms; this condition is incorporated into the solution scheme by setting 
up a flag set to zero for the attachment-line solution and to unity for * > 1. 

The rj differencing is accomplished with a combination of two schemes. For situations 
where the profile G is positive at all normal grid-point locations, the standard left-pointing 
three-point second-order scheme (which we refer to as the L scheme) is used. When 
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the profile has negative values, the “zig zag” scheme (also called the Z scheme), first 
proposed by Krause (ref. 4), is employed. For moderate crossflow situations, this scheme 
automatically satisfies the zone-of-dependence principle that is outlined by Wang (ref. 5). 
The finite-differencing formula that combines the above two schemes can be written as 


— ^1 fi,j + fafij-l + bifij—2 + bifi-ij + 65/1-1, y+l 


(131) 


The coefficients take these values for the L scheme: 

61 = (A/?j - Ar/J^/A; 62 = -Ar/^/A; 63 = 64 = 65 = 0 

A = AjfoAj/j_i(Aiy + A^_i); A rjj = rj h] - A»y_i = i/.-j- 1 - 7 7*,i-2 

The coefficients take these values for the Z scheme: 

&i = ~ mj - 1 + 

62 = -^(»?ij - ??*,>— 1); 63 =0; 64 = -^( 7 /i-ij+i - 7 /,-_ij); 65 = -64 


(132) 


(133) 


At rj = 0 (left boundary) or at rj = (right boundary), no d/dr] terms are present 
in the LISW equations; this condition is incorporated into the solution scheme by setting 
up a flag set to zero for the boundaries j = 1 and j = nylim. The flag is set to unity in 
the fully three-dimensional region. As stated earlier in section “Quasi-Two-Dimensional 
Equations for Initial Conditions,” a factor w = ufa) is used to blend the LISW solution to 
the fully three-dimensional solution at the boundary-adjacent points. The short notation 
{bF} that is used in equation (104) is thus obtained directly from equations (131)-(133). 
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4. INVISCID INTERFACE 


The three-dimensional boundary-layer solution procedure is based on the specifica- 
tion of the edge quantities u e , v e , and T e on a surface grid defined in the coordinate 
directions £ and ??. In addition, the computation of the edge density p e requires the spec- 
ification of the inviscid pressure P or the pressure coefficient C p (for flows that involve 
a shock between the free stream and the attachment line). In the general case, for a 
nonorthogonal boundary-layer grid, the metric quantities hi, h- 2 , and g u are also assumed 
to be given. Further, the edge value of viscosity p e is computed from T e with the Suther- 
land formula (equation (33)). The above quantities are referenced to the free-stream 
quantities (£/£,, P^, and T^) and the reference length Z 4 . 

If we assume that a negligible interaction occurs between the viscous and inviscid 
regions, the edge conditions can be obtained by solving the three-dimensional Euler 
equations on a sufficiently fine mesh. In some cases such as low-speed flow, a 
potential panel code may be substituted in place of the Euler solver, after which an 
interface procedure is necessary to process the inviscid results to express them in the 
form required by the boundary-layer code. Specifically, this procedure involves (a) the 
accurate location of the inviscid attachment line, (b) the generation of the boundary-layer 
grid which originates from the attachment line on the upper or lower surface, (c) the 
calculation of the edge velocities in surface grid-oriented directions, and (d) the output 
of quantities in the required form for the boundary-layer code. 

Two approaches are used to calculate the edge velocities (u e > v e) and edge tem- 
perature T e . The first approach is to interpolate the inviscid pressure distribution onto 
the boundary-layer grid and then calculate the edge velocities and temperature by solv- 
ing the limiting equations of the boundary-layer equations at the boundary-layer edge. 


30 



These limiting equations (henceforth called the BL— EDGE equations) are hyperbolic and 
can be solved with a marching method that is analogous to the boundary-layer solution 
procedure. The source term in these equations is the pressure gradient in the two di- 
rections £ and t). The second method involves the interpolation of all the required edge 
quantities from the inviscid grid to the boundary-layer grid. 

The first procedure is, in principle, more consistent with the boundary-layer solution 
method; however, the solution of the BL-EDGE equations may be slightly different from 
the solution of the Euler equations. This mismatch in the edge quantities from the two 
solutions is attributed to (a) the terms dropped in the BL— EDGE equations, based on the 
boundary-layer assumption, and (b) variations in the finite-differencing schemes. In one 
case, we have the edge velocities computed from the Euler equations, whereas, in the 
other case, the edge velocities are computed from the limiting boundary-layer equations, 
based on the Euler pressure gradient. The BL-EDGE equations also assume that the 
edge total enthalpy is a constant, which is not necessarily true for the Euler solution. 
Note also that accurate enforcement of the condition dPjdx = 0 at the attachment line 
is difficult from a coarse inviscid grid, which may necessitate that the interpolation near 
the attachment line be linear rather than spline to ensure a negative pressure gradient 
at the attachment line. The interface program includes an option to calculate the edge 
conditions by either method. 

4.1 Attachment-Line Relocation 

After the surface pressure distribution is obtained from the inviscid calculation, the 
initial location of the attachment line is obtained by scanning for a maximum pressure in 
the vicinity of the leading edge of the wing. However, note that the true attachment-point 
location may be located within a bandwidth of one grid point on either side. The true 
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location of the attachment point is where the surface velocity in the direction normal to 
the attachment line is equal to zero. 


In this procedure, the Cartesian velocity components of the Euler solution are con- 
verted to surface-oriented velocities that correspond to a boundary-layer surface grid 
generated from the initial attachment-line location. This velocity conversion is based on 
equations (1 3)— (1 7). In general, the component u e at the initial attachment point will be 
nonzero. This point is then relocated in the positive or negative direction, depending on 
the value and sign of u e and the local estimated value of du e jds\. For example, for the 
upper surface at an angle of attack, a relocation in the positive direction means that the 
point is moved to include more of the lower surface; this is done when u e has a positive 
value. After the points are relocated with this logic, a new boundary-layer surface grid 
is generated from these attachment-line locations. New values of velocity components 
are obtained by interpolation. This procedure is repeated until the value of u e at each at- 
tachment point is less than a specified tolerance value. A parameter u ; a < is used to relax 
the relocation displacement and to ensure that the iterated locations remain within the 
grid-point bandwidth mentioned above. Upon convergence, the pressure is interpolated 
with a spline routine from the inviscid grid to the final boundary-layer grid. 

4.2 Edge Values by Interpolation 

Edge velocities and temperature are obtained by spline interpolation. For a fine 
inviscid grid, this method usually produces smooth edge conditions comparable to the 
solution from the BL-EDGE equations. 

4.3 Edge Values from BL-EDGE Equations 

The three-dimensional boundary-layer equations (eqs. (19) and (20)), when applied 
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at the edge of the boundary layer, result in the BL-EDGE equations and are given as 
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Further, H e is assumed to be constant. With P specified from the inviscid solution, the 
following equation provides closure for this hyperbolic set of equations: 


T e = H e 1 + 
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Pe = P/T e 

The solution to the above system is obtained with a discretization that is identical to the 
full three-dimensional boundary-layer equations in the two surface directions. With the 
abbreviated notations for the £ and 7 directions, a discretization with Newton linearization 
yields the system 

^ (2a\u e + {au e } + b\Ve^j ^{b\u e + {ki e }) 

+ {<*v e }) (yO,\u f + 2biVe + { bv c } 





(137) 


where 


ri - W!j:{~^ Px + ^tr p y] ~ + c 25 u e v e + c 26 vl) 

r 2 = ^fcr{%?f Px ” Py } ~ ( Cs4U e + C 35U e V e + C 36 vl) 

The terms marked by an underline apply to the three-dimensional region only. These 
terms are set to zero per the LISW condition at j = 1 and j = ny. Note that the P y terms 
are retained for LISW, however. 

The solution is obtained by inversion of the system, which is followed by an iterative 
update for nonlinear terms. 
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5. BL3D EXAMPLE CASES 


Two test cases are presented here. The first test case is that of a moderately swept 
(A = 33°) tapered wing in subsonic flow. The streamwise cross sections of this wing 
correspond to NACA 0012; however, for generality, calculations of metrics and other 
parameters are done with the assumption that the wing is defined in terms of discrete 
coordinates. The second test case is a highly swept (A = 70°) wing in supersonic 
flow, similar in planform to that of the F16XL aircraft, with cross sections defined in 
discrete coordinates. For both cases, the inviscid results are obtained by solving the 
Euler equations. The computer code CFL3D (ref. 6) is used for this purpose, with the 
viscous terms set to zero for the Euler calculation. An interface routine processes the 
results and feeds the resulting edge conditions to BL3D. For validation (case 1 only), 
the results obtained from BL3D are compared with the results from the thin-layer Navier- 
Stokes code, which is also obtained with CFL3D. For the supersonic wing case, we 
present comparisons of the BL3D results with the solution from a conical swept wing 
boundary-layer code developed by Kaups and Cebeci (ref. 3). In addition, runs are 
made with a uniform suction distribution (case 1 only). These results are also compared 
with the corresponding Navier-Stokes solution. 

5.1 Geometry and Conditions for Case 1 

The planform of this wing is a trapezoid for which the root chord is 1 ft and the 
leading-edge sweep is constant at 32.73°. The wing has a span of 2 ft and the trailing- 
edge sweep is constant at 18.88°, which results in a tip chord of 0.398 ft. The tip-section 
leading edge is at x'* = 1.286 ft and the tip-section trailing edge is at x'* = 1.684 ft. 
The streamwise cross section of the wing corresponds to the NACA 0012 section. The 
free-stream conditions are = 0.5, o - 2°, P ( * - 2116 psf, and = 520°R. The wing 
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is assumed to be symmetrical about the root chord plane. 

5.2 Euler Solution for Case 1 

The surface distribution used in the Euler grid consists of constant percent chord lines 
and constant percent span lines. In the spanwise direction, the grid has 41 points on the 
wing surface, which corresponds to a span distance of 0.05 ft between grid lines. In the 
chordwise wraparound direction, the grid has 257 points, and in the wall-normal direction 
the grid has 49 points that are stretched exponentially. The grid is further extended into 
the wake region and off the tip of the wing. The Euler computation is done with the code 
CFL3D. The results from the calculation are obtained at the centers of the grid cells; 
thus, 40 cell centers exist in the spanwise direction. Let us denote these locations by 
the symbol j( INV). To avoid the region very close to the symmetry plane and the wingtip, 
we restrict our analysis to the region 7 < j(INV) < 35 , which corresponds to 0.325 < y'* ^ 
1 . 725 . The middle location of this region is atj(INV) = 21 , which corresponds to y'* = 1 . 025 . 

Figure 2(a) shows the surface distribution of the Euler grid-cell center points (note 
that only every fourth point is plotted in the chordwise direction for clarity). Figure 2(b) 
shows the corresponding boundary-layer grid on the upper surface. The boundary-layer 
grid originates from the attachment line, and the ^-coordinate is measured in terms of 
the surface arc length, which is normalized by a length such as the local chord length or, 
in the present case, by the maximum arc length. The spanwise coordinate y is defined 
as the local span distance from the symmetry plane. 

The inviscid results used in the interface routine are the three Cartesian components 
of the velocity on the wing surface, the inviscid wall density, and the temperature. 
The pressure coefficient on the surface can be calculated from the above. Alternately, 
the pressure coefficient can be specified, and the edge temperature can be calculated 
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(assuming that the edge density is given). Figure 3 shows the variation of the pressure 
coefficient obtained from the Euler solution at three span locations where j(INV) = 7, 21, 
and 35. Note that the effect of the taper is to create a favorable pressure gradient in 
the spanwise direction. 

5.3 Euler-BL3D Interface for Case 1 

The details of the procedure for relocating the inviscid attachment line are presented 
in Figures 4-7. Figure 4 shows a plot of the iterated attachment-line locations. This plot 
is in terms of (i / , £ ) of these surface points and corresponds to a view from upstream of 
the wing leading edge (not to scale). The initial attachment-line location that corresponds 
to the peak in the pressure coefficient near the leading edge is shown in solid symbols. 
The inviscid grid points on the upper and lower sides of this line are shown as dashed 
lines. These lines bound the uncertainty on the true attachment-line location because of 
the coarseness of the inviscid grid. Also shown in Figure 4 are the iterated locations of 
the attachment line. As explained in “Attachment-Line Relocation,” the iteration is based 
on a relocation strategy such that the local inviscid velocity vector is exactly tangential 
to the attachment line. Note that the line that joins the final attachment-line locations is 
much smoother than the original line. 

Figure 5 shows the velocity u e interpolated at the attachment-line location that 
corresponds to successive iterations. Note that the initial values of u e are relatively 
high (+0.04) and that the objective of driving u e to near zero (< 0.0001) is achieved 
in seven iterations. Figure 6 shows the corresponding variations of the velocity v e on 
the attachment line. Here again, the final v e variation is smooth. Figure 7 shows the 
resulting variation of C p on the attachment line. 

The boundary-layer surface grid is generated from the attachment line on the upper 
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or lower side. In the present case, 40 points are generated in the chordwise direction 
between the attachment line and the local 5 percent chord location. Beyond this point, 
the boundary-layer grid coincides with the inviscid grid, and no interpolation is required. 
The stretching of the surface distribution near the attachment line is done so that the grid 
blends smoothly with the inviscid grid (at the 5 percent chord location in the present case). 

After the generation of the surface grid, the inviscid pressure is interpolated to the 
boundary-layer grid. The boundary-layer calculation is restricted to the region that is 
bounded by the j(INV) = 7 and j{ INV) = 35 spanwise locations. The regions outside 
these limits are unsuitable; the locally infinite-swept wing assumption does not apply in 
these regions because of the proximity of the flow to the symmetry plane or the wingtip. 
The boundary-layer grid notation denoted by j(BL3D) is thus based on the left boundary 
of j(INV) = 7. In other words, the j(BL3D) = 1 location is identical to the j(INV) = 7 location. 

Figure 8 shows the chordwise variation of the interpolated C p values on the boundary- 
layer grid at three spanwise locations. The interpolation is accomplished with a spline 
routine. In the present case, because the inviscid grid has good resolution near the 
leading edge and the inviscid results are smooth in this region, no smoothing is required. 
However, in the absence of the above, a smooth spline interpolation may be required. 
In this event, the amount of smoothing and the resulting interpolated pressures must 
be carefully monitored. If the streamwise pressure gradient at the attachment line is 
not negative, a marching of the boundary-layer solution away from the attachment line 
may not be possible. Spline and other higher order interpolation methods are likely to 
introduce nonnegative pressure gradients at the attachment line. Hence, in some cases, 
a locally linear interpolation may be required near the attachment line to ensure negative 
pressure gradients. 

All physical distances, such as the root chord length, are normalized by a reference 
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length L l a . Further, the boundary-layer grid x- (and £-) coordinate is defined as the local 
chordwise arc length divided by the local maximum arc length up to the trailing edge. The 
boundary-layer grid y- (and also y) coordinate is defined as the local spanwise distance. 
With this definition, the metric quantities can be computed with equations (1)-(3). For 
the present case, the variation of hi, h 2 , and gn in the chordwise direction at the j(BL3D) 
= 1 location is shown in Figure 9. 

As explained in “Edge Values by Interpolation,” the edge velocities can be obtained 
by direct interpolation. In this method, the three Cartesian velocity components from 
the inviscid solution are interpolated. Subsequently, the edge velocities in the x and y 
directions are obtained by an inversion of the system given by equations (1 3)— (1 7). The 
edge temperature is also interpolated to complete the interface to the BL3D program. 

In the second method, the BL-EDGE equations are solved with the input pressure 
distribution as outlined in “Edge Values from BL-EDGE Equations.” Figure 10 shows a 
comparison of the edge velocities obtained by either method at the j(BL3D) = 1 location. 
Note that the two results agree well. Figure 1 1 shows contour plots of the edge velocity 
u e on the upper surface of the wing that are obtained from the two methods. A slight 
difference in the variation is caused by the fact that the BL-EDGE solution assumes 
a locally swept-wing condition at ;(BL3D) = 1 and j(BL3D) = 29. Another cause for 
this difference is that the finite differencing used in the Euler solution is different from 
that used in the BL-EDGE solution. However, for the present case, both results are 
acceptable. Figure 12 shows the corresponding comparison of the edge velocity v e . The 
edge temperature variations also compare well. 

5.4 BL3D Solution for Case 1 

The three-dimensional boundary-layer solution follows the sequence below: 
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(a) Input boundary- layer edge data. The input consists of the boundary-layer grid 
dimensions (nx, ny)\ the coordinates (x, y); the metric quantities hi, h 2 and gu\ and 
edge values u e> v e , and T e . Further, if the streamline between the free stream and 
the boundary-layer attachment line contains a shock (nonisentropic region), then the 
input of C p or P is required to calculate the edge density p e . Other inputs are the 
free-stream conditions, the reference length, and other parameters that pertain to the 
solution procedure. 

(b) Setup of initial profiles. Solution profiles that correspond to the similarity solution 
for a flat plate or wedge are used as initial profiles to start the locally infinite, swept 
attachment-line solution at i = 1, ; - 1 of the boundary-layer grid. 

(c) Setup of edge coefficients. This setup is based on the edge conditions and 
gradients. The coefficients A{, B„ Q, and D, are calculated. 

(d) LISW solution. This solution is obtained at the two boundaries j - 1 and j = nylim 
by marching in the x direction. The terms that involve derivatives in the y direction are 
set to zero for this case. Further, at i = 1, the attachment-line equations are solved. 

(e) Solution of the three-dimensional region. Given the initial plane solution and 
the side boundary solution, the three-dimensional region can now be solved with the L 
scheme or the Z scheme for the span derivatives. A switching from the L scheme to the 
Z scheme occurs when the profile of G has a negative element. 

Depending on the pressure distribution, the solution region can be restricted to (1, 
nxlim), (1, nylim). After the convergence of the solution at each ( i , j) point, quantities 
such as the boundary-layer thickness, the skin-friction coefficient, and the crossflow 
Reynolds number are calculated. 
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5.5 BL3D Results for Case 1 


Comparison of the BL3D solution is made with profiles that are obtained from the 
solution of the NS equations with the code CFL3D. This code was run on a grid with the 
same surface distribution as the Euler grid, but with 81 points in the wall-normal direction. 
The grid stretching was designed to include about 30 to 40 points in the boundary layer 
of the flow. The flow was assumed to have an abrupt transition to turbulence at the 25 
percent location on the wing. This assumption ensures that the flow remains attached 
and thereby avoids the numerical problems caused by the laminar separated regions. 
The profile comparisons are restricted to locations upstream of the 25 percent chord 
station. Because of the large amount of data, the comparison plots are presented for a 
few locations in the flow that are representative of the entire flow solution. 

Figures 13-18 compare of the BL3D and NS profiles at the j(BL3D) = 1 plane. 
We present comparisons at six locations in the chordwise direction. These locations 
approximately correspond to chord locations of 1, 2, 3, 6, 12, and 23 percent. Shown 
are the profiles of spanwise velocity v, chordwise velocity u, and temperature T. At 
these locations, the spanwise velocity profile assumes different shapes with inflection and 
reversal regions. The NS profiles are shown in open square symbols. The BL3D profiles 
at this j location are the solution of the LISW equations. In spite of this assumption, 
very good agreement is obtained until the 23 percent chord location. Note that the edge 
velocities used as input for the BL3D computation are the interpolated values from the 
inviscid code. The edge velocities that are computed from the NS solution are slightly 
different than those from the Euler solution. This difference is the main contributor to 
the lack of better agreement between the two solutions. Also, some differences are 
attributable to the fact that the profiles are not compared at exactly the same location. 
The temperature profiles at the 12 percent chord location and beyond are different 
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presumably because the boundary-layer interaction becomes more significant in this 
region. Furthermore, some differences near the attachment point are caused by the fact 
that the attachment point from the NS solution is shifted slightly compared with the Euler 
location. Overall, the agreement is satisfactory and validates the BL3D results. 

The program also calculates the crossflow within the boundary layer relative to the 
local edge streamline direction. In the present report, crossflow is defined as negative 
when pointed toward the wing root. In Figure 19, we compare the streamwise velocity 
profiles u s at six representative locations with the corresponding NS profiles. Figure 20 
shows the corresponding crossflow velocity profiles of v s . The crossflow that is predicted 
by the BL3D code is slightly larger than the NS solution in the negative crossflow region. 
The crossflow Reynolds number RecF is defined as 


Recp = 


; 5,max 
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where u* max is the maximum absolute value of the crossflow velocity v* and <$q.i is the 
normal distance at which the v* profile decays to less than 10 percent of u* max (when 
scanned from the edge down to the wall). This parameter has a strong correlation to 
the growth of crossflow instabilities in a three-dimensional boundary layer. Figure 21 
shows the variation of this parameter in the chord direction at the j(BL3D) = 1 location 
compared with the NS solution. 

Figures 22-30 present the profiles at the j(BL3D) = 15 location at the ?/ = 1.025 
spanwise station. The overall comparison is good, although in some locations differences 
in profiles exist mainly because of the fact that the edge conditions from the NS and 
Euler results are different. 

Figure 31 shows the contours of the boundary-layer thickness and the skin-friction 
coefficient in the x direction obtained from the BL3D calculation. The contours are smooth 
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and blend smoothly with the infinite swept-wing solutions at the j boundary locations. 
Note that the flow is close to laminar separation as indicated by the near-zero values 
of C ftX at £ = 0.25. 

Figure 32 shows the contours of the crossflow Reynolds number and the maximum 
absolute percent crossflow on the upper surface. The values of Re C F for this case are 
under 100, and the maximum crossflow reaches a maximum of about 12 percent. 

5.6 Results With Suction for Case 1 

Solutions were obtained from the NS and BL3D solvers with boundary-layer suction. 
A constant amount of suction q a (equal to 0.0005) was assumed. Figure 33 shows a 
comparison of the resulting solution profiles at the y(BL3D) = 15 plane. The profiles with 
no suction are also shown for comparison. Figure 34 shows the resulting Re C F values; 
a substantial reduction is produced in the crossflow with suction. 

The attachment-line Reynolds number Reg is defined as 
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where 6* m is the momentum thickness. This parameter is important because of 
attachment-line stability considerations. Figure 35 shows a comparison of the Reg values 
at the attachment line both with and without suction. The comparison with the values 
obtained from the NS solution is satisfactory. 

5.7 Geometry and Conditions for Case 2 

Case 2 is the boundary-layer flow on a supersonic wing of 70° sweep at a free- 
stream Mach number of 1.6 and at an angle of attack of 0°. The planform is similar 
to the wing of the F16XL aircraft. The other input free-stream conditions correspond 
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to an altitude of 40,000 ft ( /” = 393.13 psf, = 390°R). The tree-stream Reynolds 
number is 3.06x1 0 6 per ft. 

Figure 36 shows a top view of the Euler grid used in this case, with an inset showing 
a chordwise section. The wing is assumed to be symmetric about the (x'*,z'*) plane at 
the span station j(INV) = 1, which is at a distance of 27 in. from the fuselage axis. The 
flow region of interest corresponds to the y'* range of 72.8 to 132.2 in. 

5.8 Euler Solution and BL3D Interface for Case 2 

The inviscid pressure distribution on the wing upper surface, obtained from the Euler 
solution, is shown in Figure 37. Here, the streamwise distances are shown in terms 
of x >* _ x '* Et which corresponds to the chordwise distance from the local leading-edge 
location. The variation of pressure coefficient and the wing cross section at the span 
location j{ INV) = 16 (90.2 in. from fuselage axis) is shown in Figure 38. Attached laminar 
flow does not exist beyond a x 1 * — x'£ E distance of 3 ft; hence, the calculations reported 
here are for an x'* - x'l E of less than 3 ft. 

The interface program is essentially in the same form as case 1, except for minor 
changes in reading in and manipulation of the Euler solution input. The boundary-layer 
grid is specified as containing 40 points clustered within the 2 percent chordwise location. 
The edge velocities and temperature are calculated either by direct interpolation or by 
solution of the BL-EDGE equations. Good comparisons of the edge values u e ,v ei and 
T e from the two methods were obtained, as shown in Figures 39-41 . Figure 42 shows a 
comparison at the span station ;(BL3D) = 5. In the following section, the boundary-layer 
solution at this span station will be presented in detail. 

5.9 BL3D Solution for Case 2 

Following the interface run, the present code BL3D was run in a region bounded by 
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1 < j(BL3D) < 9. The results reported here correspond to the solution at j(BL3D) = 5. 
For comparison, the Kaups-Cebeci code was also run at this section. 

Figure 43 shows the comparison of the u s ,v s ,T profiles at two chordwise locations 
close to the attachment line at the span station j(BL3D) = 5. The agreement is good, 
except for a very small reduction in the magnitude of the crossflow profile obtained 
from BL3D. Figure 44 shows the comparison of the profiles at two locations (1 ft and 
2.26 ft in surface arc length away from the attachment line). Here, the streamwise and 
temperature profiles agree well. The crossflow profile from BL3D shows slightly reduced 
crossflow. Note that because the flow becomes increasingly three dimensional away 
from the attachment line the two codes are expected to differ in solutions. 

Figure 45 shows a comparison of the resulting crossflow Reynolds number values 
from the two computations. Again, the reduced crossflow predicted by BL3D can be 
noted. 
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Figure 1 . Three-dimensional boundary-layer coordinates for a swept wing. 
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Figure 2. Cartesian grid (a) and boundary-layer grid (b) for case 1 

(In the jc direction, only every fourth grid line is plotted for clarity). 



Figure 3. Pressure coefficient variation from Euler solution for case 1 
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Figure 4. Iterative attachment-line relocation for case 1 . 
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Figure 5. Convergence of velocity u e on attachment line for case 1 . 
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Figure 6. Convergence of velocity v e at attachment line for case 1 . 
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Figure 7. Convergence of pressure coefficient at attachment line for case 1 . 
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Figure 1 3. Comparison of NS and BL3D profiles on wing upper surface at £, = 0.0068, r| = 0.0 (Case 1 ). 
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Figure 1 4. Comparison of NS and BL3D profiles on wing upper surface at £, = 0.01 41 , r| = 0.0 (Case 1 ). 
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Figure 1 5. Comparison of NS and BL3D profiles on wing upper surface at £, = 0.0248, -g = 0.0 (Case 1 ). 
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Figure 1 6. Comparison of NS and BL3D profiles on wing upper surface at % = 0.054, r| = 0.0 (Case 1 ). 
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Figure 17. Comparison of NS and BL3D profiles on wing upper surface at ^ = 0.1 21 2, ti = 0.0 (Case 1). 



Streamwise section at ;'(BL3D) = 1 plane 
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Figure 1 8. Comparison of NS and BL3D profiles on wing upper surface at % = 0.21 86, r\ = 0.0 (Case 1 ). 
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Figure 1 9. Comparison of streamwise flow profiles from NS and BL3D at j (BL3D) = 


BL3D 
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Figure 20. Comparison of crossflow profiles from NS and BL3D at j (BL3D) = 



Streamwise section at ;'(BL3D) = 15 plane 
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Figure 22. Comparison of NS and BL3D profiles on wing upper surface at £ = 0.0067, ri = 0.7 (Case 1 ). 



Streamwise section at ; (BL3D) = 1 5 plane 
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Figure 23. Comparison of NS and BL3D profiles on wing upper surface at % = 0.01 1 7, r| = 0.7 (Case 1 ). 



Streamwise section at j (BL3D) = 1 5 plane 



Figure 24. Comparison of NS and BL3D profiles on wing upper surface at = 0.0239, t| = 0.7 (Case 1 ). 



Streamwise section at j (BL3D) = 1 5 plane 
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Figure 25. Comparison of NS and BL3D profiles on wing upper surface at £ = 0.0554, i\ = 0.7 (Case 1 ). 



Streamwise section at ;(BL3D) = 1 5 plane 
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Figure 26. Comparison of NS and BL3D profiles on wing upper surface at % = 0.1224, r| = 0.7 (Case 1 ). 


Streamwise section at y(BL3D) = 15 plane 
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Figure 27 . Comparison of NS and BL3D profiles on wing upper surface at £ = 0.2097, r\ = 0.7 (Case 1 ). 
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Figure 28. Comparison of streamwise flow profiles from NS and BL3D at j (BL3D) = 1 5. 


BL3D 
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Figure 29. Comparison of crossflow profiles from NS and BL3D at j (BL3D) = 15. 
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Figure 32. Contours of Re CF and maximum percent crossflow for case 1 , upper surface. 


Streamwise section at ;(BL3D) = 15 plane 
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Figure 33. Comparison of NS and BL3D profiles with suction at £, = 0.2097, t) = 0.7 (Case 1 ). 




Section at y(INV) = 15 
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Figure 36. Euler grid for supersonic geometry, case 2. 
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Figure 37. Contours of press 
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Euler interpolation 
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Figure 39. Comparison of contours of u e from Euler and BL-EDGE solutions. 



Euler interpolation 



Figure 40. Comparison of contours of v e from Euler and BL-EDGE solutions. 
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Figure 41 . Comparison of contours of j 
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Figure 42. Comparison of BL- EDGE solution to Euler solution for case 2, upper surface. 
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Figure 44. Comparison of solution profiles for case 2, s } *= 1.0 ft and s,* = 2.26ft. 
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Figure 45. Comparison of crossflow Reynolds numbers at j( BL3D) = 5 for case 2. 



Appendix A 


Coefficients of the Linearized System of Compressible Three-Dimensional 
Boundary-Layer Equations Discretized With a Fourth-Order Pade Formula 

(1) Coefficients of Continuity Equation: 

a n — ~ i 

a l2 = c lAr( e ll )fc — i 
a 13 = — c lfc( e 12)*_i 

«i 4 = 0 

a 15 = ~ c 2k{en)k-l 

a 16 = — c 2/:( e 12)jt_i 

a 17 ~ 0 
*11 = 1 


b n = -cu(eii)* 
b n = ~cik(eu)k 
b \4 = 0 


b l5 = c 2fc( e ll)fc 
b l6 = c 2k( e l2)k 
b l7 = 0 

r l = -(?ll)jfc + (gil)*_i + + C lk {qi 2 )k -1 - C 2 k(qu)k + C 2 A;(9l3)jk_i 


q ii = w 

q\2 = A\F^ + A 2 F + A$G V -f A±G 
qis = A\L £ + A 2 L + A^Aljf + A\M 
eil = a\A\ + A 2 
ei 2 = b iA$ -f A\ 

c\k = ^ 
c 2* = 

A 0 t = C it - Cfc -1 
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(2) Coefficients of £ Momentum Equation: 


a 2i - Fk - 1 

«22 = w k - 1 “ ^lJfc ( e 21 )*— 1 - c 2fc( e 22)fc-l 

«23 = — cut(e 2 3)fc— 1 - C2jt(e 2 4)jfc— 1 

a 24 = ~ c \kB 6 

0-25 = ~h-\ ~ C2Jt(e2l)jt_l 

a 26 = ~ c 2k( e 2i)k-\ 

a 27 = ~ c 2kB& 

b k 2l = -F k 

b\ 2 - - c l/fc( e 2l)fc + C2jfc( e 22)fc 

f>23 = C 1 A; ( e 23 ) fc + C2fc(e24)fc 
^24 = ~ c lkBe 
b 2 5 = h + c 2Jk(^2l)fc 

626 = c 2*( e 23)* 

627 = c 2kB§ 

r 2 = — (^21 )* + (<72l)fc_i + 6l*(<?22)* + C\k{<l22)k-l ~ c 2k{<l2z)k + c 2t(?23)fc_i 
921 = (IL - wF ) 

<722 = Fi (F 2 )^ + B 2 (FG) v + B 3 F 2 + B\FG + B$G 2 + BqH 

<723 = 2B\{FL)^ + B 2 (FM + GL\ + 2B3FL + B^(FM + GL ) + 2 B$GM + B$I 

e 2i = ^2aj5i + 2F 3 ) F + (^b\B 2 + B^jG 

^22 = ^2ajFi + 2 Bi^jL + ^61^2 + B^jM 
e 23 = (& 1 5 2 + j B 4 )f + 25 5 G 
e 2 4 = f B 2 + Bi'j L + 2F5A/ 
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(3) Coefficients of % Momentum Equation: 


a 31 — Gk-\ 

a 32 = — c lfc( e 3l)fc — 1 ~ c 2k( e 32)k-\ 

a 33 = w k-l — Cu-(e33)fc_i — c 2 *(e 34 ) fc _i 

«34 = ~ c lkC 6 

a 35 = — c 2k{e3l)k-l 

a 36 = ~k-l ~ C2k{C3i)k-l 

a 37 ~ ~ c 2 kG& 

& = ~G t 

^32 = - c lfc(e3l)* + c 2k{^32^k 
f*33 — ~ w k ~ C 1 Jt ( e 33 ) jt + C 2A: ( e 34)it 
^34 = c lkG$ 

^35 = c 2jfc(e3l)fc 
^36 = h + c 2it( e 33)jt 
^37 = c 2 kC(> 

r 3 = ~(?3l)jt + (93l)fc_i + CiJt(?32)jb + CUt(?32)it-l ~ c 2k(<l33)k + c 2*(?33)*_i 

9si = (IM - WG ) 

932 = Ci(FG)t + C 2 (G 2 ) v + CzF 2 + C a FG + C 5 G 2 + C 6 H 

933 = C\{FM + GX) { + 2C2(GM) v + 2C 3 FZ, + C A (FM + GL) + 2C&GM + C 6 I 

€31 — 2 C 3 F + + C^j G 

e 32 =2 C 3 L+ (aiCi+C^M 

e 33 = (aiCi + C 4 ) F + (26^2 + 2C 5 ) G 

e 3 4 = (aiCi + Ca)l+ (26^2 + 2C 5 )m 
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(4) Coefficients of Energy Equation: 


041 — Hk— 1 

0 42 = — c lJt( e 4l)fc-l - c 2Ar ( e 42 )/r — 1 
a 43 = c lJt( e 43)fc — 1 “ c 2k( e i4)k-l 

044 = U’fc-l ~ c l*( e 45)*_i ~ C2/t(e46)fc_i 

a 45 = — c 2A:( e 4l)fc_i 

046 = “ c 2fc( e 43)fc_i 

0 47 = — c 2 fc ( e 45 ) fc — 1 

6 {i = -H k 

b\ 2 = — c lfc( e 4l)fc + C2k(e42)k 

643 = C 1 Jfc ( e 43 ) A: + C 2 jt(e 4 4)* 

644 = —Wk — Cufc(e45) fc + C 2 *(e 4 6 )jt 
f>45 = c 2fc( e 4l)fc 

f»46 = c 2Jfc(e43)jt 

647 — (fp)fc "f” c 2fc( e 45)jt 

r* = -(^4l)fc + (94l)/t-l + c lfc(?42)jt + c lk(942)fc_l ~ c 2*(943)* + C 2 *(943)fc_i 
941 = Up 1 ~ wH ) 

g42 = £^77)^ + D 2 {GH) v + D^FH + D 4 GH + Ds 

943 = Z>i(FJ + HL\ + 7) 2 (G7 + HM\ + Ds(FI + HL) + D^GI + HM) + D' b 

641 — { a \D\ + 7 ) 3)77 

642 = (oi7)i + 7 ) 3)7 

643 = (7>i 7 ) 2 + 7) 4 )77 

e 44 = (6 i7 ) 2 + 7 ) 4)7 

645 = (ai 7 )i + 7)3 )F + {b\D 2 + 7) 4 )G 

e 4 6 = («i7)i + Di)L + (7>i Z>2 + D\)M 
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(5) Coefficients of F-L-L Equation: 

051 — — C2k~, £it-l 

4- 1 

052 = -1 ~ C2k~, (e5l)jt_i 

4-1 

053 = -C2k~, (C52)jb_i 

4—1 

054 = — c 2k~, Bq 

4-1 


055 — “ 

056 = 0 

057 = 0 


c \k ° 2k l k ^ "*" U, )fc-l 


^51 — c 2kJ^Lk 
f>52 = 1 + C2kJ~( e Sl)k 
hi = C 2 fc^(e 52 ) fc 
f>54 = C2fe — Be 


1 


&55 = —Ci k + C2fc — ( — l' + w) t 

b 56 = 0 
657 = 0 


r 5 = - (95l)jt + 1 + Clfc(952)*. + c ljt(952)jb-l ~ c 2<;(953),t + c 2 *(953)jt_ 


951 = ^ 

952 = 4 


g 53 = y|^(F 2 ) ? + B 2 (^G% + B,F 2 + B*FG + B b G 2 + B%H 

— LI * + tv L + A\FF^ + 

e^i = ^ 2a\B\ + 2B$ + a\A\^ F A (b\B2 + B^j G + A\F je -f A^G^ 

C 52 = ( 'b 1 B 2 + B i + b 1 A^F + 2B 5 G 
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(6) Coefficients of G-M-A/' Equation: 


«61 — — C2jt 7 — Mk - 1 
* Jt-i 

062 = -C2k~, ( e 61 )fc — 1 

4-1 

«63 = -1 - C 2 k~, (e62)fc_i 

4-i 

064 — —C2k~, C*6 

4-1 

065 = 0 

066 = “Clfc - c 2k + “Ojb-1 

067 = 0 


4;i = c 2k~r^h 
4 

662 = C M~ k ( e 6l)fc 

663 = 1 + C2kJ^{ e §2)k 
f>64 = c 2kl~Cfj 

b 65 = 0 k 

f>66 = ~C\k + c 2kJ~ { — l' + w )k 
ber = 0 


r 6 = ~(<l6l)k + (%l)k-\ + c lk{%2)k + c \k(%2)k-l ~ c 2k(<l63)k + ^2* (^63 ) fc — 1 


961 = G 

962 = M 

963 = y|ci(FG)£ + C 2 (G 2 ) ji + C 3 F 2 + C±FG + C b G 2 + C 6 H 

- Ml' + wM + AiGFt + A 3 GG v j 

C6i = 2 C b F + {^a\C\ + C4 + a\A\^jG 

e62 — (oiCi + C*4^ F + ( 2 b\C 2 + 2C5 + &i^3^ G A\F^ -\- A b Gjj 
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(7) Coefficients of H-l-l Equation: 


1 


«71 = ~ c 2k~ 


-Wk - 1 


(^p)fc-l 

072 = ~ c 2kTT-, ( e 7l)jt_i 

Wfc-1 

073 = ~C2k— ( e 72)jt_l 




074 = — 1 — C 2k 

075 = 0 

076 = 0 


077 = ~C\k ~ C 2 k 


1 


(‘P/Jfc-l 


( e 73)*-l 


('p)jfc- 1 


( -/ p + ^*-1 


*71 = C 2 fc- 


1 




’(Wt 

*72 = C 2 jfcTTT-( e 7l)jt 

VP>k 

*73 = o 2 fc 77 - r (e 7 2)it 


(g 


P7jt 


*74 = 1 + C 2k 

*75 = 0 


1 

(Wi 


'( e 73)A 


*76 = 0 

*77 = — C 1 fc + c 2k , i (~^p + O') A; 

Wit 

p 7 = -(<77i)a + (97l)*_i + c lk (q 72 ) k + cifc(9 72 ) fc _j - c 2k (q 73 ) k + c 2k (q 73 ) k 
971 = H 


q 72 

973 



+ D 2 (GH\ + D 3 FH + D a GH + D 5 
-Ifr + wI + AiHFt + AiHGA 


e 7i = fajDi + i^3 + 

672 = f*i D 2 + D 4 + biA^jH 

673 = + P 3 )f+ (*iZ) 2 + ^ 4 )g + A 1 ^ + ^ 3 ^ 
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Appendix B 


User’s Manual for the BL3D Program 

The BL3D program package consists of a collection of program files, input files, and 
output files that relate to the solution of three-dimensional boundary-layer flow on a swept 
wing. In addition to the BL3D program files, the program package includes files related 
to the Euler and Navier-Stokes calculations for the two example cases and the interface 
programs. Some of these files are not directly related to the core program BL3D, but 
were developed for the validation of the code. These include Euler and Navier-Stokes 
grid and solution files as well as programs used in processing the Navier-Stokes results 
for comparison with the BL3D results. All the files have been archived on the NASA 
Langley Masstor system and can be made available per individual request. 

Table 1 at the end of this section gives a list of important computer variable names 
used in the programs and the corresponding translation to the variables used in the 
formulation as listed in “Nomenclature”. This list may be useful for customizing the 
program for a particular application. 

Program Structure 

The files are divided into subdirectories as follows. Each subdirectory contains a 
readme file with a short description of the contents of the subdirectory and instructions 
for running the programs. 

bl3d/euler/casel . Contains files that correspond to generation of the Euler 
solution for case 1 (subsonic swept wing) with the code CFL3D. The program CFL3D 
is not included. 

bl3d/interface/casei. Contains interface programs and input files used in 
processing the Euler solution for case 1 and generation of the data file required in the 
BL3D run. 
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bi3d/bl /easel. Contains the BL3D program files used for case 1. 
bl3d/ns/casel. Contains files that correspond to the Navier-Stokes solution with 
CFL3D for case 1. 

Similarly, bl3d/euler/case2, bl3d/interface/case2, bl3d/bl/case2, 
bl3d/ns/case2 are subdirectories that correspond to case 2 runs for a supersonic 
swept wing. The subdirectory bl3d/ns/case2 contains files that relate to the applica- 
tion of the Kaups-Cebeci code for case 2. 

The files in the subdirectories bl3d/euler and bl3d/ns are not of direct interest 
here. The readme files in these subdirectories and comment statements in the program 
files may be consulted for more detailed information. Presented below are details on the 
interface and boundary-layer program files for the two test cases. 

Interface Program Inputs 

The inviscid results are assumed to be available on a grid that is generally (though not 
strictly) oriented along the constant span and constant chord directions. The boundary- 
layer grid stretching, the interpolation, and the BL-EDGE solution procedure are based 
on this assumption. The surface grid is also, in general, nonorthogonal. The inviscid 
Euler results are assumed to be available in either of the two following formats: 

(a) Two separate files, one of which contains the grid points and the other, which 
contains the Euler solution at the cell center points. The actual read statement (e.g., 
from bl3d/interface/casel/rinvis . f) is as follows: 

The grid file: 

open(unit=l, file='grid' , form= ' unformatted ' , status= ' old ' ) 
read(l) jdim, kdim, idim 

read (1) ( ( (xg ( j , k, i) , j=l, jdim) , k=l , kdim) , i=l , idim) , 

( ( (yg ( j , k, i) , j=l, jdim) , k=l , kdim) , i=l, idim) , 

{ ( (zg ( j , k, i) , j=l, jdim) , k=l, kdim) , i=l, idim) 

The solution file: 
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dimension titlw(20) 

open (unit=l , file= 'rest .bin' , form= 'unformatted' , status='old' ) 
read(2) titlw, xmachw, j t, kt, it , alphw, ntr, ntime 

read (2) ( ( ( (q ( j , k, i , 1 ) , j=l, jdim-1) , k=l,kdim-l) , i=l, idim-1) , 1=1, 5) 


The above file format corresponds to output from the program CFL3D. The i , j , k 
indices in the read statement are defined differently compared to the BL3D terminology. 
Here, i refers to the span direction, and j refers to the chord wraparound direction that 
starts from the wing lower surface trailing edge or the downstream plane of the lower 
surface wake. Note that in this case the solution points are one less than the grid points 
because they are obtained at the cell centers. Following these read statements, the 
i, j indices are switched to correspond to BL3D notation and are normalized with the 
reference length refien. The coordinates and solution on the wing surface alone are 
then extracted. In some cases, the sign of the coordinate y and velocity component vc 
may have to be switched because in BL3D the coordinate system assumes a right wing. 
The variable q contains the solution vector in standard CFL3D notation (see the readme 
file for more detailed information). 

(b) A provision is included in the interface program to input the inviscid data on the 
surface only from a single ASCII file in the following format: 


open (unit=21 , f ile= ' edge . dat ' , status^ 'unknown' ) 
read (21 , * ) iv, jv 

read (21, *) ( (xc (i, j) , i=l, iv) , j=l, jv) 
read (21, *) ( (yc(i, j) , i=l, iv) , j=l, jv) 
read (21, *) ( (zc (i, j ) , i=l, iv) , j=l, jv) 
read(21, *) ( (uc (i, j) , i = l, iv) , j = 1 , jv) 
read (21, *) ( (vc(i,j) , i=l, iv) , j=l, jv) 
read (21, *) ( (wc (i , j ) , i=l , iv) , j=l, jv) 
read (21, *) ( (tc (i, j) , i=l, iv) , j=l, jv) 
read (21, *) ( (pc (i, j) , i=l, iv) , j=l, jv) 
read (21, *) ( (xs (i, j) , i=l, iv+1) , j=l, jv+1) 
read (21 , * ) ( (ys(i, j) ,i=l,iv+l) , j=l, jv+1) 
read (21, *) ( (zs (i, j ) , i=l, iv+1) , j=l, jv+1) 


The variables iv, jv refer to the surface grid cell center dimensions in the chord 
wraparound direction and the span direction. The variables xc,yc, zc are coordinates 
of the cell centers (assumed to be nondimensionalized with refien) and uc , vc , wc are 
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the Cartesian velocity components (normalized by free-stream velocity) from the inviscid 
solution at the cell centers. The variables tc,pc are nondimensional temperature and 
pressure coefficient values at the cell centers. The variables xs,ys, zs are the actual 
normalized grid locations (not used in the program except for plot output). 

The input file for the interface run is in the following format: 

b!3d/ inter face/ easel /inp.dat 


CASE 1, 

INTERFACE 

RUN 

istart 

amach 

gam 

isurf 

0.5 

1.40 

1 

0 

itel 

i te2 

j tel 

j te2 

33 

289 

1 

43 

jl 

j 2 

ref len 


7 

35 

1.0 


nep 

40 

itat 

xcp 

0.05 

uerlim 

omatt 


50 

0.0002 

0.8 

iout 

nxw 

j lw 

j2w 

100 

7 

35 

1 


Note that each input line is preceded by a line that contains a description of the 
input variables (read as a character variable and ignored). The description of the input 
variables is as follows: 


amach 

gam 

isurf 

istart 


itel, ite2, jtel, jte2 


jl/ j2 

ref len 


Mach number 

ratio of specific heats 

upper or lower surface (1 or 0) 

restart option 

=0 full run, input from 

inviscid files 'grid', 'rest.bin' 

=1 run with data read from 'edge.dat' 

=2 read inviscid files, write 'edge.dat' 
and stop 

indices corresponding to wing edges 
itel , i te2 : lower and upper surface TE 
indices {inviscid data assumed to be 
indexed in the chord-wrap-around dir- 
ection from lower TE to LE to upper TE) 
jtel,jte2; span direction wing boundaries 

limit inviscid data processing to (jl->j2) 

normalizing length, usually root chord; 
if grid units are in inches for example, 



ncp , xcp 


itat , uerlim, omatt 


nxw, j lw, j2w 


iout 


reflen can be used to be consistent with 
units 

number of points in the BL grid between 
attachment point and local chord fraction 
=xcp (30,0.05 for example for 30 points 
from attach, point to 5% chord location) 

max number or attachment point relocation 
iterations and error tolerance on ue 
(20,0.001 recommended) 

omatt is a relaxation factor for attach- 
ment point relocation iteration, <1 ; 
smaller values give slower convergence, 
but less oscillation 

BL3D output option 

output file write will be from 

i=l to i=nxw, j = j lw to j=j2w 

nxw . LE . nx ; j lw . GE . j 1 ; j 2w . LE . j 2 

if nxw is set to 0, nxw taken to be-nx 

=1, output interpolated ue,ve,te values 
=2, output ue,ve,te from 

solution of BL-EDGE equations 


Note that the istart=l option corresponds to the input of inviscid data from a single 
ASCII file as described previously. If the calculation is for the lower surface (isurf=o) 
and if istart=0, then the program reverses the i indexing direction (in effect, the lower 
surface becomes the upper surface of an inverted wing). 


Interface Program Output 

The output for the BL3D code is written to a binary file called bi .bin. The format 
for this write is as follows (see also the file blout . f ) : 


nyw= j 2 w- j lw+ 1 

write ( 10 ) nxw, nyw . 

write ( 10 ) ( (x(i, j) ,i=l,nxw) , j=:lw, j 2 w 
write ( 10 ) { (y(i, j) ,i=l,nxw) , j=Dlw, 32 w) 
write ( 10 ) ( (hl(i, j) ,i = l,nxw) , j = jlw,;j 2 w 
write ( 10 ) ( (h 2 (i, j) ,i=l,nxw) , j=jlw,j 2 w) 
write ( 10 ) ( (gl 2 (i, j) ,i=l,nxw) , j = 3 lw, ;| 2 w) 
if (iout. eq. 1 ) then . 

write ( 10 ) ( (ue(i,n) ,i=l,nxw) , 3 =d lw, d 2 w) 
write ( 10 ) ( (ve(i, 3 ) ,i = l,nxw) , j lw, ;} 2 w) 

0 1 SG 

write ( 10 ) ( (used, j) ,i = l,nxw) , j=jlw, j2w) 
write ( 10 ) ( (vse(i, j) ,i=l,nxw) , j=jlw, 32 w) 
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endi f 

write (10) ( (pb ( i , j ) , i = l , nxw) , j=jlw, j2w) 
if ( iout . eg . 1 ) then 

write (10) ( ( tb(i , j ) , i=l, nxw) , j= jlw, j2w) 

0l S0 

write (10) ( (tse(i,j) ,i=l,nxw) , j=jlw, j2w) 
endif 

write ( 10 ) ( (xb(i, j) ,i=l,nxw) , j=jlw, j2w) 
write (10) ( (yb ( i , j ) , i=l,nxw) , j=jlw, j2w) 
write (10) ( (zb(i, j ) , i=l,nxw) , j=jlw, j2w) 
write ( 10) ( (ub(i, j) ,i=l,nxw) , j=jlw, j2w) 
write (10) ( (vb(i, j) , i=l,nxw) , j=jlw, j2w) 
write ( 10 ) ( (wb ( i , j ) , i=l, nxw) , j= j lw, j2w) 
write ( 10 ) (xsave(nx, j ) , j=jlw, j2w) 

The output of xsave corresponds to the normalizing arc length values used in each 
span station j=jlw, j2w. Note that the dimensional distance from the attachment line 
in the x direction (i.e., surface arc length) is given by x(i, j) *xsave( j) *reflen. 
Other output files contain information about attachment-line iteration convergence, the 
accuracy of the interpolation, and the BL— EDGE solution. The files readme, main, f, 
as well as comment statements throughout the various subroutines, may be consulted 
for more information on these outputs. 


Running the Interface Program 


The program routines are contained in the following files: 


main.f 
attach. f 
blout . f 
inputs . f 
metrix. f 
pack. f 
ploti . f 
rinvis . f 
sdist . f 
solv3 . f 


surf . f 


In addition, the include files com, compack, the makefile makefile, and the input 
file inp.dat are also provided. A brief description of some of these files follows: 
inputs, f. Reads the input file inp.dat. 
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rinvis . f . Reads the inviscid data; extracts required data needed for the upper or 
lower surface; normalization of coordinates and flow quantities. 

attach, f. Locates the initial attachment line (for first pass only); generates a 
boundary- layer surface grid that originates from the currently defined attachment line; 
calls sdist . f for a suitable grid stretching. 

metrix.f. Calculates metric quantities hi,h2,gi2 based on currently defined 
boundary-layer grid. 

surf. f. Solves the BL-EDGE equations if this input option is selected. 

blout . f. Outputs the file bl .bin required for BL3D run. 

pack. f. Library of cubic tension spline interpolation routines. 

The include file com has the following parameter statements (which correspond to 
case 1 in this example): 

parameter ( jdim=321, kdim=49 , idim=53 ) 
parameter (ivm=2 57, jvm=43) 
parameter (nxm=3 00 , nym= jvm) 

The dimensions depend on the inviscid grid size and the boundary- layer grid size; 
see comment statements in the file com for selection of these values. 

The include file compack is used in conjunction with the spline routines in pack. f. 

The program is run as follows: choose appropriate dimensions in the parameter 
statements in the file com; remove all . o files if com is modified; then type make followed 
by a. out to run the program. The output files fort . 3 , fort .20, fort . 30, fort . 31, 
fort . 32 , fort . 50 , fort . 51 , fort . 99 are useful for plotting or checking the results. 
Please refer to the various subroutines for a description of these outputs. 

An appropriate selection of input quantities and inviscid data of good resolution are 
required for successful completion of the interface run. Run times on the Cray are only a 
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few CPU seconds so that runs can usually be made interactively. Some error messages 
and suggested fixes are given below. 


error from s/r SDIST; ipass=... 

no solution for ak found in the range . . . to . . . ' , 
change ak or dak values in s/r ATTACH or 

change input values of ncp,xcp 

The possible problem here is the use of inappropriate boundary-layer grid stretching; 
try a different combination of ncp,xcp. The inviscid grid may also be too coarse near 
the leading edge. 


Convergence to uerlim not achieved 

This problem is due to the lack of convergence in the attachment- line relocation 
loop; choose a lower value of omatt or increase itat, uerlim values. If various 
combinations do not work, inviscid data is too coarse or the solution near the attachment 
line is not smooth enough or the inviscid attachment line is too curved. 


no convergence; s/r surf 


i, 3 =• • • ... 

ii , use (ii, j ) , vse ( ii , 3 ) 


hse (ii , j ) ,pb{ii , j ) 


This problem is due to lack of convergence in the BL-EDGE equations, which is 
invariably caused by a small nonnegative dP/dx at the attachment line. The nonnegative 
dP/dx is caused by an amplification in the oscillations in pressure near the attachment 
line from the spline interpolation caused by a coarse inviscid grid. One possible fix is to 
drop back to linear interpolation for P for a few points near the attachment line. Some 
comment statements can be found in attach, f, which can be “uncommented” to drop 
back to linear interpolation. 
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Boundary-Layer Program Inputs 


The inputs to the boundary-layer program are composed of two files: 

(a) The boundary-layer edge data file bi .bin. This binary file is created by running 
the interface program and is assumed to be available in the current directory (bi /easel, 
for example). The format for this file has been described previously. This file is read 
in the subroutine rinvis. 


(b) The boundary-layer input options file inp . dat. This file is read in the subroutine 
inputs and is assumed to be in the following format (note that each input line is 
preceded by a line that contains a description of the input variables, which is read as 
a character variable and ignored): 
b!3d/bl/ easel /inp. dat 


CASE 1 : BL INPUTS 


AMACH 

PFS 

TFS 

PRL 

IWALL 

REFLEN 

IUNIT 

0.5 

2000 

520 

0.72 

0 

1.0 

0 

ZMAX 

AK 

NXLIM 

NYLIM 




8.0 

1.05 

68 

29 




ITMAX 

EPSF 

EPSG 

EPSH 

IORD 




20 l.e-6 l.e-6 l.e-06 5 


The description of the input variables is given below: 


amach 


free-stream Mach number 


pf s 
tfs 
prl 
iwall 


ref len 


free-stream static pressure (lb/ft2) 

free-stream static temperature (deg R) 

laminar Prandtl number 

wall boundary condition 
=0 adiabatic wall 

=1 wall temperature, deg R specified 
(read in 'main.f') 

=2 wall heat flux specified ,1b/ (ft. sec) 

(read in main.f) 

note: suction rates are also specified in 'main.f' 
reference length in ft. 

(similar to input in interface program, but in ft.) 
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iunit 

zmax 

ak 

nxlim 

nylim 

itmax 

epsf 

epsh 

iord 


US or SI units, =0 for US units 
SI units not currently implemented 

maximum value of transformed variable 'zeta' 
corresponding to boundary-layer edge 

stretching coefficient by which points are 
distributed away from the wall 

number of stations to march in x direction 

number of span-wise stations in y direction 

maximum number of iterations at each (i,j) point 

convergence criterion for f-prime at wall 

convergence criterion for h-prime at wall 

value of i for which scheme becomes fully second- 

order in the stream-wise direction 


Boundary-Layer Program Outputs 

The output from the boundary-layer program is output to several fort, files. The 
convergence information can be found in fort. 2 . Most of the results are output from 
subroutine phys and can be modified as required. The output quantities include 
z*, u, v, T profiles; boundary-layer thickness, displacement thickness, and skin-friction 
coefficient; streamwise and crossflow profiles; and crossflow Reynolds number. Please 
refer to subroutine phys for more details. Note that these “writes” are for each (/, 
j) location (with j as the inner loop). 

Output Format for Stability Analysis 

The output of profiles and their derivatives to first and second order, as well as 
edge-normalizing quantities, are made in the format described below. 


open (uni t=7 , f ile= ' fort . 7 ' , form= 'unformatted' , status= 'unknown' ) 
write ( 7 ) title 

Title in character*80 . 
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write (7) nxlim, nylim, reflen 


nxlim is number of chordwise stations, 
nylim is number of spanwise stations, 
reflen is reference length, dimensional. 

write (7) pf s, tf s , ufs , amach, rof s 
Quantities 7^, tf*,, 

For each profile i=l, nxlim; j = l, nylim: 

write (7 ) i,j,nz,x(i,j) ,y (i, j ) , xbs,ybs, zbs 

write (7 ) hi (i, j ) , h2 (i, j ) , gl2 (i , j ) , sis 

write (7 ) ues , ves , tes , roes , pes , res , ales, prl , del tax, the tax, rcr 
i is streamwise index, 
j is spanwise index, 
nz is no. of points in profile (constant), 
x ( i , j ) is boundary-layer coordinate jc. 
y ( i , j ) is boundary-layer coordinate y. 
xbs , ybs , zbs are Cartesian coordinates x'*, /*, z!*. 
hi ( i , j ) , h2 (i , j ) , gl2 { i , j ) are metrics h h h 2 , gn- 
sis is s\*, dimensional arc-length in the r direction, 
ues, ves, tes, roes, pes are edge quantities u e *, v e *, T e *, p e *, P e *. 

For i=l, du*Jds\ is output instead of u e *. 
res is local Reynolds number based on u e *, v e * . 
ales is reference length by which normal distance is scaled, 
ales = 1* = for i=1 - 
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ales = l* = y/^fs\Ju* for i>l. 
prl is laminar Prandtl number. 

del tax is boundary-layer thickness (99 percent) based on u, dimensional, 
thetax is momentum thickness based on u, dimensional, 
rcr is radius of curvature, dimensional. 


Profile outputs: 


write(7) (zp(k) ,k=l,nz) 

zp is z p = (z*/l*)\ note l* definition for i=l and i>l. 


This box only for i>l: 

write (7 ) (up (k) , k=l , nz) 
write (7) (upz (k) , k=l,nz) 

write (7) (upzz(k) , k=l,nz) 

up is u*/u e *\ upz is du p /dz p ] upzz is d 2 u v /dz 2 


write(7) (vp(k) ,k=l,nz) 
write (7 ) (vpz (k) , k=l , nz ) 

write (7 ) (vpzz (k) , k=l , nz ) 

vp is v*/v e * if i=l, v*/u e * if i > 1; 

vpz, vpzz are derivatives defined similar to upz, upzz. 

write(7) (tp(k) ,k=l,nz) 
write (7) (tpz(k) ,k=l,nz) 

write{7) (tpzz(k) ,k=l,nz) 

tp is T*/T*\ tpz , tpzz are defined similar to upz, upzz. 
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write(7) (wp(k) ,k=l,nz) 
write ( 7 ) (wpz (k) , k=l , nz ) 

write (7) (wpzz (k) , k=l,nz) 

wp is w*/v* e if i=l, w*/u* for i>l; wpz, wpzz are defined similar to upz,upzz. 


This box only for i>2; output of streamwise derivatives: 


write(7) (duds(k) ,k=l,nz) 
write (7) (dvds(k) ,k=l,nz) 
write(7) (dtds(k) ,k=l,nz) 
write (7) (dwds(k) ,k=l,nz) 

duds=(du*/ds\)/{U^/r) 


dvds =(dv*/dsl)/(UUn 

dtds =(dT*/ds\)/(U*jn 
dwds =(dw*/ds\)/{U^/l*) 


These output writes are implemented in subroutine out in the subdirectory, 
bl3d/bl/casel (and also bl3d/bl/case2). 

Running the Boundary-Layer Program 

The program routines are contained in the following files: 


blk. f 
der. f 
etac . f 
main . f 
phys . f 
rinvis . f 
sumk . f 
cl3phi . f 
derk . f 
initial . f 
nextep . f 
read7 . f 
solve . f 
updte . f 


B-13 




conv. f 
edge . f 
inputs . f 
out . f 
ref . f 
sura, f 

In addition, the include files com and compack, the makefile makefile, and the 
input file inp . dat are also provided. A brief description of important subroutines follows: 

inputs, f. Reads input file inp. dat. 

ref . f. Calculates reference quantities. 

rinvis. f. Reads boundary-layer edge data bl .bin. 

initial . f. Sets up initial profile at the attachment point i = 1 , j = 1 . 

main. f. Calls various routines; specification of wall temperature or wall heat flux; 
specifies wall suction rate. 

edge. f. Specifies edge coefficients depending on boundary point or interior point; 
calls cl3phi . f , etac . f to calculate some additional edge parameters. 

etac . f . Sets up the y differencing coefficients based on the local crossflow direction 
(L scheme or Z scheme). 

blk.f. Sets up and solves the (7x7xnz) block tridiagonal system. 

updte. f. Updates nonlinear quantities 6, 6', l, /', D 5 , and / at k = 1 for the current 
iteration. 

updte . f . Checks convergence. 

phys.f. Outputs dimensional quantities, boundary-layer properties; calculates 
crossflow. 

nexstep. f: storing of variables prior to stepping in /. 

out. f: output of quantities and profiles for boundary-layer stability analysis. 

The include file com has the following parameter statement (for easel in this 
example): 
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parameter (nx=100 , ny=29 , nz=61) 

The nx and ny values are as written out from the interface program. The value of 
nz corresponds to the number of points in the boundary-layer normal direction. 

The program is run with input files inp.dat and bi.bin. Appropriate dimensions 
in parameter statement in the file com are selected; then remove all .o files if com 
was modified. Type make, followed by a. out to run the program. Runs can be made 
interactively for moderate grids. The program execution stops at the (i,j) location where 
the boundary layer separates; this results from the calculation of negative temperature 
in subroutine update. An error message is printed out, and the incipient separation 
condition can be checked by the local skin-friction coefficient. 


Table 1 


List of Important Program Variables 


Computer variable Symbol 


a, b, c 

a l,m' 7 l,m 

al , a2 , a3 , a4 

Mi M, M, M 

adl , ad2 , ad3 

a 1, d2, «3 

ak 

k s 

al 

l 

alp 

i 

alpr 

I/Pr 

alprp 

I'/Pr 

amach 

Moo 

amue 

Pe 

amuel 

Pe at (i,j) 

amuf s 

* 

ftoc 

amur 

fl r 

bl , b2 , b3 , b4 , b5 , b6 

5i, 5 2 , 53 , 54 

55, 56 

b3h, b4h 

53, 5 4 


Description 

elements of linearized system (eq. (120)) 

coefficients of transformed continuity 
equations (eqs. (47)-(50)) 

( direction coefficients (eq. (129)) 

stretching constant for ( distribution (eq. 
(128)) 

ratio C pn)/(pe^e ) 
dl/d c 

Mach number 

edge viscosity, nondimensional, (eq. (33)) 
local value of n e 

free-stream viscosity, dimensional 

reference viscosity function, 
nondimensional (eq. (33)) 

coefficients (eqs. (55)- (60)) 
coefficients (eq. (94)) 
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Computer variable 

Symbol 

Description 

b3 1 , b4 1 , b5 1 , b6 1 

B3, B\ , B$ , Be 

coefficients (eq. (91)) 

bdl , bd2 , bd3 , bd4 , 
bd5 

bi, bi, 63. 64, 65 

?/ direction differencing coefficients (eq 
(132)) 

bet 

7 /( 7 - 1 ) 


bet2 

(7 — l )/2 


01,02,03,04,05, c6 

Cu C 2 , 03,0-4, 

C$, Ce 

coefficients (eqs. (62) -(67)) 

c4h, c5h 

CaA 

coefficients (eq. (95)) 

c3t , c4t, c5t , c6t 

&AAA 

coefficients (eq. (92)) 

cl3 

Ci 3 

metric coefficient (eq. ( 6 )) 

cl32 

Ci 3 2 


cl3phix 

( C\z<f>lh\) x 

metric coefficient (eq. (48)) 

cl3phiy 

{Cn<j>v r /(h2U e )} 7 y 

, metric coefficient (eq. (50)) 

c24, c25 , c26 

C24) C25, C26 

coefficients (eq. (19)) 

c2 6x 

dC 2 e/dx 

coefficient (eq. (80)) 

c34, c35, c36 

Cm, C35, C36 

coefficients (eq. ( 20 )) 

cbetal 

cos (/3) IJ 

angle between grid lines (eq. (3)) 

cdl, cd2 

1 1 

2’ 12 

constants (eq. ( 101 )) 

chil, chi2 , chi3 

Xl,X 2 >X 3 

coefficients (eq. (87)) 
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Computer variable Symbol 


dl,d2,d3,d4,d5 

Dl,D-2,Dz,D4, 

D 5 

d5p 

dD^/dC 

d3 t , d4t 

7*3, -C*4 

delfp, delgp, 
delhp 

dS 5 ,dS 6 ,dS 7 

delfpmax, 

delgpmax, 

delhpmax. 

max. of 

dS 5 ,dS 6 ,dS 7 

del tax 


duds 

(du e /dx) l=l 

el , e2 , e3 , e4 

Ei, E'2, E 3 , E 4 

epsf , epsg, epsh 


f 

F 

fl, f2 


fe 

{du e jdx) i=l 

fel 


fey 

{d 2 u e /dxdy) i= 

fp 

F' 

fpl , fp2 



Description 

coefficients (eqs. (68)-(69)) 

coefficients (eq. (96)) 

solution vector elements (eq. (108)) 

convergence indicators 

boundary-layer thickness based on F 
reaching 99 percent value, dimensional 

attachment-line velocity gradient 

coefficients (eq. (87)) 

convergence criterion for change in 
f,g',ti at wall {k= 1) 

u/u e 

F at i - 1 and i - 2 
attachment-line velocity gradient 
local value of fe 
used in eq. (76) 

d(u/u e )/d( 

F* at i - 1 and / - 2 
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Computer variable 


Symbol 


g 

fl, f 2 

gi, g2 

gl2,gl21 

gi2x 

gl2y 

gam 

gel 

gex, gey 

gp 

gpl, gp2 
h 
hi 
hll 

hlx, hly 

h2 

h21 

h2x, h2y 
he 


Description 


G 


8n 


dgnjdx 

dgn/dy 

7 


G' 


H 


hi 


h , 

^■ 1 , 1 ? ^ 1 ,y 

hi 


hi 

h2,x, ^2,y 


Hk=ke 


v/v T 

F at i - 1 and / - 2 
G at / - 1 and / - 2 
metric coefficient 


local value of Gk=ke 

x ,y gradients of G k=ke 

d{v/v r )ld(i 

G' at i — 1 and i - 2 

nondimensional total enthalpy (eq. (24)) 

metric coefficient 

local value of metric coefficient 

*,>’ gradients of h\ 

metric coefficient 

local value of metric coefficient 

x,y gradients of h 2 

boundary-layer edge value of H 
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Computer variable Symbol 


Description 


hel 


local value of H 

hfs 

H'oo 

free-stream value of total enthalpy (eq. 
(23)) 

hhl , hh2 


H at i - 1 and i - 2 

hp 

H' 

dH/d( 

hpl , hp2 


H 1 at i - 1 and i - 2 

hse 


H e from BL-EDGE equations 

i 

i 

streamwise index 

iat 


attachment-line iteration count 

idim 


Euler grid maximum dimension 

ile 


i index of leading edge 

iord 


value of i at which scheme is fully second 
order in x 

iout 


flag to output edge values by interpolation 
or from BL-EDGE solution 

istart 


flag for interface run mode 

isurf 


flag for upper or lower surface 

it 

superscript n 

boundary-layer iteration count 

itat 


maximum number of attachment-line 
iterations 

itel , i te2 


wing trailing edge indices 
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Computer variable 

itmax 

iunit 

iv 

ivm 

iwall 

D 

jl, j2, j3 , j4 , j5 
j lw, j2w 
jdim 

j tel, j te2 
jv 

jvm 

kdim 

lflagx 

If lagy 


Symbol Description 


maximum number of boundary-layer 
iterations 

flag for U.S. units or S.l. units 

number of points in streamwise direction, 
from interface run 

maximum dimension for iv 

flag for wall boundary condition 

spanwise index 

neighboring points of (i , j ) used in y 
differencing 

spanwise indices that correspond to 
interface output 

Euler grid maximum dimension 

wing symmetry plane and tip indices 

number of points in spanwise direction 
from interface run 

maximum dimension for jv 

Euler grid maximum dimension 

flag to signify attachment-point solution or 
three-dimensional solution 

flag to signify LISW or general 
three-dimensional case 
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Computer variable Symbol 


Description 


ncp 


nx 

nxlim 

nxm 

nxw 

ny 

nylim 

nym 

nz 

nzl 

omatt 

omega 

pb 

pc 

pfs 

phil 


C P 

C P 

4 > 


number of points in boundary-layer grid 
between attachment point and point at 
which local chord fraction = xcp 

streamwise dimension of boundary-layer 
grid 

value of i<nx limiting march in x 

maximum dimension for nx 

number of streamwise stations written out 
from interface code 

spanwise dimension of boundary-layer grid 
value of j (<ny) limiting sweep in y 
maximum dimension for ny 
number of points in normal direction 
nz-l 

relaxation factor u a for attachment-line 
iteration 

blending factor between LISW and 
three-dimensional solutions 

pressure coefficient on boundary-layer grid 

pressure coefficient on Euler grid 

free-stream static pressure, dimensional 

local value of <p 
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Computer variable 

Symbol 

Description 

pi 

7T 


prl 

Pr 

laminar Prandtl number 

q 

q 

absolute velocity, nondimensional 

qel 


local value of absolute edge velocity 

ref len 

L'oo 

reference length 

ref Is 

Reoo 

free-stream Reynolds number based on 
L*oo (eq. (29)) ' 

refs 


free-stream Reynolds number per unit 
length 

roe 

Pe 

boundary-layer edge density, 
nondimensional 

roel 

Pe 

local value of boundary-layer edge density 

rof s 

Poo 

free-stream density, dimensional 

rstar 

R* 

gas constant 

sll 

Sl 

nondimensional arc length 

tb 

T e 

nondimensional boundary-layer edge 
temperature 

tc 


nondimensional temperature from Euler 
solution 

te, tel 

T e 

nondimensional boundary-layer edge 
temperature 

tfs 

rp* 

1 oo 

free-stream dimensional temperature 
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Computer variable Symbol 


Description 


thel , the2 

M' 

temperature ratio and gradient (eqs. (54) 
and (89)) 

the tax 


momentum thickness based on F, 
dimensional 

til, ti2 


temperature profiles at i — 1 and i - 2 

tr 

T r 

reference temperature (eq. (33)) 

trs 

T r * 

reference temperature, dimensional (eq. 
(32)) 

tse 


edge temperature from the solution of the 
BL-EDGE equations 

tst 


streamwise derivative of T 

tw 

Tw 

wall temperature, nondimensional (eq. 
(119)) 

ub 


Cartesian velocity component, 
nondimensional 

uc 


Cartesian velocity component from Euler 
solution, nondimensional 

ue 

u e 

streamwise edge velocity 

uel 

U e 

local value of streamwise edge velocity 

uerlim 


error tolerance on streamwise velocity u e 
at attachment line 

uex 

du e /dx 

streamwise velocity gradient 

uey 

dujdy 

spanwise velocity gradient 
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Computer variable Symbol 


Description 


uf s 

U'oc 

freestream velocity, dimensional 

uil , ui2 


velocity profiles at i - 1 and i - 2 

url 


= 0 on attachment line, = 1 , away from 
attachment line; used to compute q 

use 


edge velocity from the solution of the 
BL-EDGE equations 

ustr 


velocity profile in the edge streamline 
direction 

vb 


Cartesian velocity component, 
nondimensional 

VC 


Cartesian velocity component from Euler 
solution, nondimensional 

vcross 


crossflow velocity (orthogonal to edge 
streamline) 

ve 

Ve 

edge velocity in span direction, 
nondimensional 

vel 

Ve 

local value of v e 

vil , vi2 


velocity profiles at i - 1 and i - 2 

vrl 

V r 

reference velocity used in the definition of 
G (eq. (43)) 

vrx , vry 


derivatives in x,y of v r 

vse 


edge velocity from the solution of the 
BL-EDGE equations 

w 

W 

transformed normal velocity (eq. (44)) 
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Computer variable Symbol 


Description 


wb 

wc 

wil,wi2 

X/Y 

xb, yb, zb 

xc , yc , zc 
xcp 
xnorm 

z 

zcon 
zil , zi2 
zphys 

zphysx 

zphysy 


Cartesian velocity component, 
nondimensional 

Cartesian velocity component from Euler 
solution, nondimensional 

velocity profiles at i — 1 and / - 2 

.x,y streamwise boundary-layer coordinates 

x' Cartesian coordinate of boundary-layer 

grid 

Cartesian coordinate of Euler grid 

percent chord value at which i=ncp 

normalizing arc length, nondimensional 

transformed boundary-layer normal 
coordinate 

coefficient to calculate zphys from z 
zphys at i — 1 and i — 2 

z* boundary-layer grid normal coordinate, 

dimensional 

streamwise gradient of z* 
spanwise gradient of z* 
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